Meta-analysis of avian responses and butterfly eyespots

1 Setting ups

1.1 Load packages and custum function

Load packages and custum function. Custom function is used for calculating effect size (lnRR) and its variance from raw data.

Code

pacman::p_load(ape,
               GoodmanKruskal,
               DT,
               dtplyr,
               ggokabeito,
               ggcorrplot,
               here, 
               knitr,
               tidyverse,
               patchwork,
               metafor,
               miWQS,
               multcomp,
               orchaRd,
               parallel)

# function for calculating effect size (lnRR)
effect_lnRR <- function(dt) {
  # replace 0 in "C_mean", "T_sd", "C_sd", "C_proportion" with each minimum values
  # if proportion too extreme value, replace minimum value (only "T_proportion")

  dt <- dt %>%
  mutate(
    C_mean = ifelse(C_mean == 0, 0.04, C_mean),
    T_sd = ifelse(T_sd == 0, 0.01, T_sd),
    C_sd = ifelse(C_sd == 0, 0.05, C_sd),
    C_proportion = ifelse(C_proportion == 0, 0.01, C_proportion),
    T_proportion = ifelse(T_proportion < 0.01, 0.01, T_proportion),
    T_proportion = ifelse(T_proportion == 1, 0.9, T_proportion)
  )
    
  # copy dataset for adding effect size and its variation (lnRR /lnRR_var) column
  dt1 <- dt %>%
    mutate(
      lnRR = NA,
      lnRR_var = NA
    )

  for (i in seq_len(nrow(dt1))) {
    Tn <- dt1$Tn[i]
    Cn <- dt1$Cn[i]
    T_mean <- dt1$T_mean[i]
    C_mean <- dt1$C_mean[i]
    T_proportion <- dt1$T_proportion[i]
    C_proportion <- dt1$C_proportion[i]
    T_sd <- dt1$T_sd[i]
    C_sd <- dt1$C_sd[i]
    Response <- dt1$Response[i]
    Measure <- dt1$Measure[i]
    Study_design <- dt1$Study_design[i]
    Direction <- dt1$Direction[i]

    # continuous data - using escalc() (metafor package)
    if (Response == "continuous" & Study_design == "independent") {

      effect <- escalc(
        measure = "ROM",
        n1i = Tn, n2i = Cn,
        m1i = T_mean, m2 = C_mean,
        sd1i = T_sd, sd2i = C_sd,
        var.names = c("lnRR", "lnRR_var")
      )

      dt1$lnRR[i] <- effect$lnRR
      dt1$lnRR_var[i] <- effect$lnRR_var
    } else if (Response == "continuous" & Study_design == "dependent") {

      effect <- escalc(
        measure = "ROMC",
        ni = (Tn + Cn) / 2,
        m1i = T_mean, m2 = C_mean,
        sd1i = T_sd, sd2i = C_sd,
        ri = 0.5,
        var.names = c("lnRR", "lnRR_var")
      )

      dt1$lnRR[i] <- effect$lnRR
      dt1$lnRR_var[i] <- effect$lnRR_var
    }

    # proportion data (no sd values)
    else if (Response == "proportion1" & Study_design == "independent") {
      
    T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )

      # transform proportion value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro1 <- log(T_proportion / C_proportion)
      lnRR_var_pro1 <- (1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn) +
        1 / (C_proportion^2 * Cn))

      dt1$lnRR[i] <- lnRR_pro1
      dt1$lnRR_var[i] <- lnRR_var_pro1

    } else if (Response == "proportion1" & Study_design == "dependent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )

      # transform proportion value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro1 <- log(T_proportion / C_proportion)
      lnRR_var_pro1 <- (1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn) +
        1 / (C_proportion^2 * Cn)) -
        2 * 0.5 * sqrt((1 / sqrt(8))^2 * (1 / (T_proportion^2 * Tn))) *
          sqrt((1 / sqrt(8))^2 * (1 / (C_proportion^2 * Cn)))

      dt1$lnRR[i] <- lnRR_pro1
      dt1$lnRR_var[i] <- lnRR_var_pro1
    }

    # proportion data (exist sd values)
    else if (Response == "proportion2" & Study_design == "independent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )
      
      # transform proportion mean value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_SD <- sqrt(T_sd^2 / (4 * (T_proportion) * (1 - (T_proportion))))
      C_SD <- sqrt(C_sd^2 / (4 * (C_proportion) * (1 - (C_proportion))))

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro2 <- log(T_proportion / C_proportion)
      lnRR_var_pro2 <- (T_SD)^2 * (1 / (T_proportion^2 * Tn)) +
        (C_SD)^2 * (1 / (C_proportion^2 * Cn))

      dt1$lnRR[i] <- lnRR_pro2
      dt1$lnRR_var[i] <- lnRR_var_pro2

    } else if (Response == "proportion2" & Study_design == "dependent") {
      
      T_proportion <- replace(
        T_proportion, Direction == "Decrease",
        (1 - T_proportion[Direction == "Decrease"])
      )
      C_proportion <- replace(
        C_proportion, Direction == "Decrease",
        (1 - C_proportion[Direction == "Decrease"])
      )
      
      # transform proportion mean value
      asin_trans <- function(proportion) {
        trans <- asin(sqrt(proportion))
        trans
      }

      T_SD <- sqrt(T_sd^2 / (4 * (T_proportion) * (1 - (T_proportion))))
      C_SD <- sqrt(C_sd^2 / (4 * (C_proportion) * (1 - (C_proportion))))

      T_proportion <- asin_trans(T_proportion)
      C_proportion <- asin_trans(C_proportion)

      # calculate lnRR and lnRR variance
      lnRR_pro2 <- log(T_proportion / C_proportion)
      lnRR_var_pro2 <- (T_SD)^2 * (1 / (T_proportion^2 * Tn)) +
        (C_SD)^2 * (1 / (C_proportion^2 * Cn)) -
        2 * 0.5 * sqrt((T_SD)^2 * (1 / (T_proportion^2 * Tn))) *
          sqrt((C_SD)^2 * (1 / (C_proportion^2 * Cn)))

      dt1$lnRR[i] <- lnRR_pro2
      dt1$lnRR_var[i] <- lnRR_var_pro2
    }
  }

  return(dt1)
}

1.2 Prepare datasets for analysis

Read datasets and phylogeny

First, we used a predator dataset to check whether we need to consider phylogeney. If so, we will use predator and prey datasets for following meta-analysis.

Code

dat_all <-  read.csv(here("data/all_21092023.csv"), header = T, fileEncoding = "CP932")
dat_pred <- filter(dat_all, Dataset == "predator")
dat_prey <- filter(dat_all, Dataset == "prey")

datatable(dat_all, caption = "all dataset", options = list(pageLength = 10, scrollX = TRUE))

Table 1 - Predator + prey datasets

Code

trees <- read.nexus(here("data/bird_phy.nex"))
plot(trees[1])

Figure 1— Phylogenetic tree of bird species

Calculate lnRR and lnRR variance and prepare datasets for analysis

Code

# predator dataset
dt_pred <- effect_lnRR(dat_pred)
dt_pred$Obs_ID <- 1:nrow(dt_pred)
dt_pred$Log_diameter <- log(dt_pred$Diameter_pattern)
dt_pred$Log_area <- log(dt_pred$Area_pattern)
dt_pred$Log_background <- log(dt_pred$Area_background)

# prey dataset
dt_prey <- effect_lnRR(dat_prey)
dt_prey$Obs_ID <- 1:nrow(dt_prey)
dt_prey$Log_diameter <- log(dt_prey$Diameter_pattern)
dt_prey$Log_area <- log(dt_prey$Area_pattern)
dt_prey$Log_background <- log(dt_prey$Area_background)

# all dataset
dt_all <- effect_lnRR(dat_all)
dt_all$Obs_ID <- 1:nrow(dt_all)
dt_all$Log_diameter <- log(dt_all$Diameter_pattern)
dt_all$Log_area <- log(dt_all$Area_pattern)
dt_all$Log_background <- log(dt_all$Area_background)

hist(dt_all$lnRR)

Code

hist(dt_all$lnRR_var)

Code

# use vcalc to calculate variance-covariance matrix
VCV_pred <- vcalc(vi = lnRR_var, 
                  cluster = Cohort_ID, 
                  subgroup = Obs_ID, 
                  data = dt_pred)

VCV_prey <- vcalc(vi = lnRR_var, 
                  cluster = Cohort_ID, 
                  subgroup = Obs_ID, 
                  data = dt_prey)

VCV <- vcalc(vi = lnRR_var, 
             cluster = Cohort_ID,
             obs = Obs_ID,
             rho = 0.5,
             data = dt_all)

dt_pred <- dt_pred %>%
  mutate_if(is.character, as.factor)
dt_prey <- dt_prey %>%
  mutate_if(is.character, as.factor)
dt_all <- dt_all %>%
  mutate_if(is.character, as.factor)

summary(dt_pred)

          Study_ID          Cohort_ID      Shared_control_ID      Year     
 Brilot_2009  :40   Jones_1980_1 :16   Jones_1980_a :16      Min.   :1957  
 Jones_1980   :16   Brilot_2009_1:10   Brilot_2009_a:10      1st Qu.:2009  
 Vallin_2011  :10   Brilot_2009_2:10   Brilot_2009_b:10      Median :2009  
 Hossie_2015  : 8   Brilot_2009_3:10   Brilot_2009_c:10      Mean   :2005  
 Olofsson_2015: 7   Brilot_2009_4:10   Brilot_2009_d:10      3rd Qu.:2013  
 Vallin_2010  : 5   Vallin_2011_1:10   Vallin_2011_a:10      Max.   :2016  
 (Other)      :31   (Other)      :51   (Other)      :51                    
    Country         Data_source Data_location              Bird_species
 Canada : 8   table. 1    :21   p.187  :24    Cyanocitta_cristata: 2   
 Finland:14   fig. 3      :15   p.186  :16    Emberiza_sulphurata: 1   
 India  : 2   fig. 2      :12   p.214  :16    Ficedula_hypoleuca : 2   
 Sweden :24   text, fig. 3:10   p.1663 :10    Gallus_gallus      :48   
 UK     :67   text        : 9   p.13   : 5    Parus_caeruleus    :19   
 US     : 2   fig. 4a     : 8   p.1419 : 4    Parus_major        : 3   
              (Other)     :42   (Other):42    Sturnus_vulgaris   :42   
        Bird_common_name            Prey_species
 blue jay       : 2      Caligo_martia    : 2   
 blue tit       :19      Inachis_io       : 6   
 common starling:42      Junonia_almana   : 3   
 domestic fowl  :48      Lasiommata_megera:11   
 great tit      : 3      none             :88   
 pied flycatcher: 2      Pieris_rapae     : 2   
 yellow bunting : 1      Saturnia_pavonia : 5   
                   Prey_common_name        Response        Measure  
 emperor moth              : 5      continuous :75   latency   :61  
 European cabbage butterfly: 2      proportion1:16   number    :15  
 none                      :88      proportion2:26   proportion:39  
 owl butterfly             : 2                       time      : 2  
 peacock butterfly         : 6                                      
 peacock pansy butterfly   : 3                                      
 wall brown butterfly      :11                                      
    Direction                                        Detail_result
 Decrease:24   latency to attack                            :10   
 Increase:81   proportion of time spent near eyespot        : 9   
 Neutral :12   latency to approach the food bowl            : 8   
               latency to first movement                    : 8   
               propotion of time spent furthest from eyespot: 8   
               proportion of time spent on foodbowl         : 7   
               (Other)                                      :67   
                                                                                                                 Note_result 
                                                                                                                       :107  
 including wing flaps and/or jumping                                                                                   :  1  
 mean time in seconds from the beginning of a trial to when a blue tit first visited the log on the floor (first visit):  3  
 the time from landing on the perch in front of the butterfly until the attack was launched                            :  4  
 total number of day 1 to 5, reversed value                                                                            :  1  
 weak and intense alarm call                                                                                           :  1  
                                                                                                                             
   Treatment_stimulus Number_pattern Diameter_pattern  Area_pattern    
 conspicuous:  7      Min.   :1.00   Min.   : 2.20    Min.   :  3.801  
 eyespot    :110      1st Qu.:1.00   1st Qu.: 4.50    1st Qu.: 15.904  
                      Median :2.00   Median :20.00    Median : 87.715  
                      Mean   :1.94   Mean   :15.29    Mean   :254.622  
                      3rd Qu.:2.00   3rd Qu.:25.00    3rd Qu.:489.697  
                      Max.   :5.00   Max.   :31.28    Max.   :490.874  
                                                                       
 Area_background Shape_pattern   Control_stimulus
 Min.   :  225   circle :  6   camouflage: 10    
 1st Qu.: 1000   eyespot:110   no pattern:107    
 Median : 3200   stripe :  1                     
 Mean   : 5851                                   
 3rd Qu.:13175                                   
 Max.   :26400                                   
                                                 
                                                                                                                                                                                                       Note_stimulus
                                                                                                                                                                                                              :22   
 single wing; painted over eyespot on about half of the butterflies (ﾔeyespot painted overﾕ) and sham-painted an equally large area just beside the eyespot on the rest of the butterflies (ﾔeyespot visibleﾕ): 6   
 ambiguous eyespot; based on a photographof owl eyes with a contrasting light iris and dark pupil; with alarm call                                                                                            : 5   
 ambiguous eyespot; based on a photographof owl eyes with a contrasting light iris and dark pupil; with sparrowhawk call                                                                                      : 5   
 ambiguous eyespot; based on a photographof owl eyes with a contrasting light iris and dark pupil; with threat call                                                                                           : 5   
 ambiguous eyespot; based on a photographof owl eyes with a contrasting light iris and dark pupil; with white noise                                                                                           : 5   
 (Other)                                                                                                                                                                                                      :69   
      Type_prey                 Shape_prey        Diet_bird     Bird_sex 
 artificial:88   abstract_butterfly  :14   grainivore  :  1   both  :50  
 real      :29   abstract_caterpillar:16   Invertebrate:  2   female:12  
                 abstract_prey       :58   omnivore    :114   male  :16  
                 butterfly           :29                      NA's  :39  
                                                                         
                                                                         
                                                                         
     Bird_age        Tn              Cn            T_mean         
 adult   :48   Min.   : 4.00   Min.   : 4.00   Min.   :   0.0004  
 chick   :38   1st Qu.: 8.00   1st Qu.: 8.00   1st Qu.:   0.7656  
 juvenile: 4   Median : 9.00   Median : 9.00   Median :  16.1000  
 NA's    :27   Mean   :14.73   Mean   :14.26   Mean   : 119.1540  
               3rd Qu.:15.00   3rd Qu.:15.00   3rd Qu.: 153.5000  
               Max.   :89.00   Max.   :89.00   Max.   :1180.4878  
                                               NA's   :16         
     C_mean             MT_mean             MC_mean          
 Min.   :   0.0400   Min.   :-556.1000   Min.   :-2115.0000  
 1st Qu.:   0.5079   1st Qu.:  -6.8415   1st Qu.:  -13.6000  
 Median :   9.3000   Median :   0.7544   Median :    0.5543  
 Mean   : 110.1728   Mean   : -43.7981   Mean   : -121.8786  
 3rd Qu.: 154.2000   3rd Qu.:   0.9522   3rd Qu.:    0.9174  
 Max.   :2115.0000   Max.   :   0.9996   Max.   :    0.9574  
 NA's   :16          NA's   :81          NA's   :81          
      T_sd               C_sd            T_proportion      C_proportion    
 Min.   :  0.0100   Min.   :   0.0500   Min.   :0.01000   Min.   :0.04262  
 1st Qu.:  0.5183   1st Qu.:   0.2859   1st Qu.:0.05258   1st Qu.:0.08262  
 Median : 16.9814   Median :   4.7434   Median :0.27152   Median :0.26477  
 Mean   : 82.2725   Mean   :  78.2657   Mean   :0.30879   Mean   :0.29854  
 3rd Qu.: 74.7486   3rd Qu.:  56.2141   3rd Qu.:0.49500   3rd Qu.:0.50000  
 Max.   :731.8482   Max.   :1052.1749   Max.   :0.82609   Max.   :0.79618  
 NA's   :16         NA's   :16          NA's   :75        NA's   :75       
      Study_design     Dataset    Mean_median T_mean_median      
 dependent  :64    predator:117   mea   : 1   Min.   :   0.0004  
 independent:53                   mean  :99   1st Qu.:   0.7656  
                                  median: 1   Median :  12.0000  
                                  NA's  :16   Mean   : 104.2470  
                                              3rd Qu.: 104.0000  
                                              Max.   :1180.4878  
                                                                 
 C_mean_median       Type_of_variance_statistic      T_se          
 Min.   :   0.0000   IQR : 1                    Min.   :  0.00000  
 1st Qu.:   0.6559   SD  : 3                    1st Qu.:  0.09462  
 Median :   9.3000   SE  :28                    Median :  4.88000  
 Mean   :  96.8310   SE  : 1                    Mean   : 25.83973  
 3rd Qu.:  61.0000   SEM :68                    3rd Qu.: 24.68328  
 Max.   :2115.0000   NA's:16                    Max.   :182.92683  
                                                NA's   :20         
      C_se          T_variance_q1  T_variance_q3 C_variance_q1 C_variance_q3
 Min.   :  0.0000   Min.   :77.5   Min.   :246   Min.   :29    Min.   :36   
 1st Qu.:  0.0903   1st Qu.:77.5   1st Qu.:246   1st Qu.:29    1st Qu.:36   
 Median :  1.5000   Median :77.5   Median :246   Median :29    Median :36   
 Mean   : 26.7360   Mean   :77.5   Mean   :246   Mean   :29    Mean   :36   
 3rd Qu.: 19.8747   3rd Qu.:77.5   3rd Qu.:246   3rd Qu.:29    3rd Qu.:36   
 Max.   :372.0000   Max.   :77.5   Max.   :246   Max.   :29    Max.   :36   
 NA's   :20         NA's   :116    NA's   :116   NA's   :116   NA's   :116  
                                                                                                                                          Note   
                                                                                                                                            :72  
  authors only reported total sample size (left eyespots and right eyespots), so the sample size was equally divided between the two stimuli: 8  
 Mann-Whitney U test; p < 0.002                                                                                                             : 4  
 ｱ 1 SEM                                                                                                                                    : 4  
 Calculated by AM                                                                                                                           : 3  
 (Other)                                                                                                                                    :24  
 NA's                                                                                                                                       : 2  
      lnRR             lnRR_var             Obs_ID     Log_diameter   
 Min.   :-2.15769   Min.   :0.0008295   Min.   :  1   Min.   :0.7885  
 1st Qu.:-0.07696   1st Qu.:0.0107871   1st Qu.: 30   1st Qu.:1.5041  
 Median : 0.08779   Median :0.0333636   Median : 59   Median :2.9957  
 Mean   : 0.27498   Mean   :0.1249595   Mean   : 59   Mean   :2.4180  
 3rd Qu.: 0.54880   3rd Qu.:0.1655345   3rd Qu.: 88   3rd Qu.:3.2189  
 Max.   : 4.16777   Max.   :1.3885984   Max.   :117   Max.   :3.4430  
                                                                      
    Log_area     Log_background  
 Min.   :1.335   Min.   : 5.416  
 1st Qu.:2.767   1st Qu.: 6.908  
 Median :4.474   Median : 8.071  
 Mean   :4.574   Mean   : 7.819  
 3rd Qu.:6.194   3rd Qu.: 9.486  
 Max.   :6.196   Max.   :10.181  
                                 

Code

summary(dt_prey)

                 Study_ID            Cohort_ID               Shared_control_ID
 Hossie_2013         :25   Hossie_2013_13 : 13   Hossie_2013_g        : 13    
 Hossie_2012         :24   Hossie_2012_19 :  6   Stevens_2007_b       :  8    
 Stevens_2007        :20   Ho_2016_13     :  2   deWert_2012_a        :  7    
 Ho_2016             :18   Ho_2016_14     :  2   Hossie_2012_j        :  6    
 Stevens&Hardman_2008:14   Ho_2016_15     :  2   Stevens_2007_d       :  6    
 Stevens&Cantor_2009 :12   Vallin_2010_2_1:  2   Stevens&Cantor_2009_a:  6    
 (Other)             :40   (Other)        :126   (Other)              :107    
      Year           Country      Data_source                Data_location
 Min.   :2003   Canada   :49   fig. 4   :32   supplementary Material:54   
 1st Qu.:2008   Finland  : 2   fig. S2  :30   p.387                 :12   
 Median :2012   Germany  : 2   fig. 3   :18   p.1223                :10   
 Mean   :2011   Singapore:18   table. S1:18   p.529                 : 9   
 3rd Qu.:2013   Sweden   : 8   fig. 5   : 8   p.1222                : 8   
 Max.   :2022   UK       :71   fig. 7   : 8   p.27                  : 7   
                US       : 3   (Other)  :39   (Other)               :53   
             Bird_species        Bird_common_name            Prey_species
 Ficedula_hypoleuca:  2   blue tit       :  8     Bicyclus_anynana :  2  
 none              :143   none           :143     Inachis_io       :  2  
 Parus_caeruleus   :  8   pied flycatcher:  2     Mycalesis_perseus: 18  
                                                  none             :116  
                                                  Saturnia_pavonia :  8  
                                                  Thyatira batis   :  7  
                                                                         
                   Prey_common_name        Response         Measure   
 Bicyclus anynana          :  2     continuous : 19   number    : 19  
 dingy bush brown butterfly: 18     proportion1:131   proportion:134  
 emperor moth              :  8     proportion2:  3                   
 none                      :116                                       
 peacock butterfly         :  2                                       
 Thyatira batis            :  7                                       
                                                                      
    Direction              Detail_result
 Decrease:  6   attacked          :  6  
 Increase:128   escaped           :  1  
 Neutral : 19   number of attacked: 19  
                survival          :127  
                                        
                                        
                                        
                                        Note_result    Treatment_stimulus
                                              :124   conspicuous:77      
 2 cm length                                  :  1   eyespot    :76      
 4 cm length                                  :  1                       
 6 cm length                                  :  1                       
 AM culculated the survival rate from raw data:  2                       
 calculated by AM; (total-predated)/total     : 24                       
                                                                         
 Number_pattern   Diameter_pattern  Area_pattern    Area_background 
 Min.   : 1.000   Min.   : 1.700   Min.   :  2.27   Min.   : 129.6  
 1st Qu.: 2.000   1st Qu.: 4.500   1st Qu.: 15.90   1st Qu.: 625.0  
 Median : 2.000   Median : 6.000   Median : 28.27   Median : 913.9  
 Mean   : 2.562   Mean   : 6.759   Mean   : 35.30   Mean   : 891.2  
 3rd Qu.: 2.000   3rd Qu.: 8.380   3rd Qu.: 43.01   3rd Qu.:1110.6  
 Max.   :11.000   Max.   :15.000   Max.   :176.71   Max.   :2784.7  
                                                                    
   Shape_pattern   Control_stimulus
 circle   :54    bh        :  1    
 diamond  : 5    bhv       :  1    
 eyespot  :76    bv        :  1    
 lectangle:14    bvv       :  1    
 others   : 4    feces     :  3    
                 no pattern:146    
                                   
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Note_stimulus
 eyespots-defensive posture vs. no eyespots-defensive posture ;models were based loosely on the late instars of Papilio canadensis and Papilio glaucus caterpillars                                                                                                                                                                                                                                                                                                                                                                                 : 12  
 single, paper models of the dingy bush brown butterfly; model which has reduced wings and eyespot diameter of approximately 6_mm. A mirror image of the butterfly was created, and a band (approx. 6_mm in length and 5_mm in width) with the same colour as the butterfly's body was added between the two images connecting the mid-section of the body. The butterfly images were then printed in colour on normal paper using a colour laser printer. The paper butterflies were cut out and dipped in paraffin wax to render them waterproof. :  9  
 solid-eyespots vs. solid-no eyespots ;models were based loosely on the late instars of Papilio canadensis and Papilio glaucus caterpillars                                                                                                                                                                                                                                                                                                                                                                                                         :  9  
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    :  8  
 countershaded-eyespots vs. countershaded-no eyespots ;models were based loosely on the late instars of Papilio canadensis and Papilio glaucus caterpillars                                                                                                                                                                                                                                                                                                                                                                                         :  8  
 models were based loosely on the late instars of Papilio canadensis and Papilio glaucus caterpillars                                                                                                                                                                                                                                                                                                                                                                                                                                               :  7  
 (Other)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            :100  
      Type_prey                  Shape_prey        Diet_bird   Bird_sex  
 artificial:116   abstract_butterfly  :64   Invertebrate:  2   NA's:153  
 real      : 37   abstract_caterpillar:52   omnivore    :  8             
                  butterfly           :37   NA's        :143             
                                                                         
                                                                         
                                                                         
                                                                         
     Bird_age         Tn                Cn              T_mean      
 nestling:  1   Min.   :  7.708   Min.   :  7.708   Min.   :0.1148  
 NA's    :152   1st Qu.: 20.000   1st Qu.: 20.000   1st Qu.:0.4656  
                Median : 60.000   Median : 60.000   Median :0.7374  
                Mean   : 63.020   Mean   : 63.026   Mean   :1.0794  
                3rd Qu.: 96.000   3rd Qu.: 96.000   3rd Qu.:1.2864  
                Max.   :160.000   Max.   :160.000   Max.   :6.1300  
                                                    NA's   :134     
     C_mean          MT_mean           MC_mean             T_sd       
 Min.   :0.1552   Min.   :-6.1300   Min.   :-4.3900   Min.   :0.1403  
 1st Qu.:0.3357   1st Qu.:-1.2048   1st Qu.:-0.9296   1st Qu.:0.6219  
 Median :0.4496   Median :-0.5582   Median :-0.3664   Median :0.9192  
 Mean   :0.8256   Mean   :-0.6946   Mean   :-0.4626   Mean   :1.2667  
 3rd Qu.:1.1137   3rd Qu.:-0.1148   3rd Qu.:-0.1552   3rd Qu.:1.8102  
 Max.   :4.3900   Max.   : 0.8822   Max.   : 0.9584   Max.   :5.4821  
 NA's   :134      NA's   :128       NA's   :128       NA's   :131     
      C_sd         T_proportion      C_proportion         Study_design
 Min.   :0.0864   Min.   :0.01598   Min.   :0.0100   dependent  :  1  
 1st Qu.:0.1422   1st Qu.:0.19510   1st Qu.:0.1100   independent:152  
 Median :0.2532   Median :0.31679   Median :0.1896                    
 Mean   :0.4709   Mean   :0.36705   Mean   :0.2783                    
 3rd Qu.:0.4555   3rd Qu.:0.49938   3rd Qu.:0.4125                    
 Max.   :3.6547   Max.   :0.90000   Max.   :0.8333                    
 NA's   :131      NA's   :19        NA's   :19                        
 Dataset    Mean_median  T_mean_median    C_mean_median   
 prey:153   mean  :  1   Min.   :0.1085   Min.   :0.1424  
            median: 18   1st Qu.:0.4496   1st Qu.:0.3137  
            NA's  :134   Median :0.6934   Median :0.4112  
                         Mean   :1.0349   Mean   :0.7883  
                         3rd Qu.:1.2320   3rd Qu.:1.0391  
                         Max.   :6.1300   Max.   :4.3900  
                         NA's   :134      NA's   :134     
 Type_of_variance_statistic      T_se            C_se        T_variance_q1    
 IQR : 18                   Min.   :0.040   Min.   :0.0200   Min.   :0.07075  
 SE  :  1                   1st Qu.:0.055   1st Qu.:0.0350   1st Qu.:0.28640  
 NA's:134                   Median :0.065   Median :0.0450   Median :0.47294  
                            Mean   :0.185   Mean   :0.1225   Mean   :0.52460  
                            3rd Qu.:0.195   3rd Qu.:0.1325   3rd Qu.:0.72800  
                            Max.   :0.570   Max.   :0.3800   Max.   :1.22170  
                            NA's   :149     NA's   :149      NA's   :135      
 T_variance_q3    C_variance_q1    C_variance_q3                    Note    
 Min.   :0.1651   Min.   :0.1136   Min.   :0.2096                     :146  
 1st Qu.:0.6584   1st Qu.:0.2594   1st Qu.:0.4316   Big asymmetry 1   :  1  
 Median :0.9294   Median :0.3159   Median :0.5600   Big asymmetry 2   :  1  
 Mean   :1.1200   Mean   :0.4762   Mean   :0.8185   Control with spots:  1  
 3rd Qu.:1.6352   3rd Qu.:0.7712   3rd Qu.:1.3267   Mid asymmetry 1   :  1  
 Max.   :2.3216   Max.   :0.9481   Max.   :1.6880   Mid asymmetry 2   :  1  
 NA's   :135      NA's   :135      NA's   :135      (Other)           :  2  
      lnRR             lnRR_var            Obs_ID     Log_diameter   
 Min.   :-0.81549   Min.   :0.002201   Min.   :  1   Min.   :0.5306  
 1st Qu.:-0.01671   1st Qu.:0.011480   1st Qu.: 39   1st Qu.:1.5041  
 Median : 0.21564   Median :0.022422   Median : 77   Median :1.7918  
 Mean   : 0.24895   Mean   :0.058391   Mean   : 77   Mean   :1.8249  
 3rd Qu.: 0.46508   3rd Qu.:0.047242   3rd Qu.:115   3rd Qu.:2.1258  
 Max.   : 2.24770   Max.   :0.434523   Max.   :153   Max.   :2.7081  
                                                                     
    Log_area      Log_background 
 Min.   :0.8197   Min.   :4.864  
 1st Qu.:2.7666   1st Qu.:6.438  
 Median :3.3420   Median :6.818  
 Mean   :3.2987   Mean   :6.673  
 3rd Qu.:3.7614   3rd Qu.:7.013  
 Max.   :5.1745   Max.   :7.932  
                                 

Code

summary(dt_all)

         Study_ID            Cohort_ID       Shared_control_ID      Year     
 Brilot_2009 : 40   Jones_1980_1  : 16   Jones_1980_a : 16     Min.   :1957  
 Hossie_2013 : 25   Hossie_2013_13: 13   Hossie_2013_g: 13     1st Qu.:2009  
 Hossie_2012 : 24   Brilot_2009_1 : 10   Brilot_2009_a: 10     Median :2010  
 Stevens_2007: 20   Brilot_2009_2 : 10   Brilot_2009_b: 10     Mean   :2009  
 Ho_2016     : 18   Brilot_2009_3 : 10   Brilot_2009_c: 10     3rd Qu.:2013  
 Jones_1980  : 16   Brilot_2009_4 : 10   Brilot_2009_d: 10     Max.   :2022  
 (Other)     :127   (Other)       :201   (Other)      :201                   
      Country       Data_source                 Data_location
 UK       :138   fig. 3   : 33   supplementary Material: 54  
 Canada   : 57   fig. 4   : 32   p.187                 : 24  
 Sweden   : 32   fig. S2  : 30   p.186                 : 16  
 Singapore: 18   table. 1 : 23   p.214                 : 16  
 Finland  : 16   fig. 2   : 19   p.387                 : 12  
 US       :  5   table. S1: 18   p.1223                : 10  
 (Other)  :  4   (Other)  :115   (Other)               :138  
             Bird_species        Bird_common_name            Prey_species
 none              :143   none           :143     none             :204  
 Gallus_gallus     : 48   domestic fowl  : 48     Mycalesis_perseus: 18  
 Sturnus_vulgaris  : 42   common starling: 42     Saturnia_pavonia : 13  
 Parus_caeruleus   : 27   blue tit       : 27     Lasiommata_megera: 11  
 Ficedula_hypoleuca:  4   pied flycatcher:  4     Inachis_io       :  8  
 Parus_major       :  3   great tit      :  3     Thyatira batis   :  7  
 (Other)           :  3   (Other)        :  3     (Other)          :  9  
                   Prey_common_name        Response         Measure   
 none                      :204     continuous : 94   latency   : 61  
 dingy bush brown butterfly: 18     proportion1:147   number    : 34  
 emperor moth              : 13     proportion2: 29   proportion:173  
 wall brown butterfly      : 11                       time      :  2  
 peacock butterfly         :  8                                       
 Thyatira batis            :  7                                       
 (Other)                   :  9                                       
    Direction                                 Detail_result
 Decrease: 30   survival                             :127  
 Increase:209   number of attacked                   : 19  
 Neutral : 31   latency to attack                    : 10  
                proportion of time spent near eyespot:  9  
                latency to approach the food bowl    :  8  
                latency to first movement            :  8  
                (Other)                              : 89  
                                                                                                                 Note_result 
                                                                                                                       :231  
 calculated by AM; (total-predated)/total                                                                              : 24  
 the time from landing on the perch in front of the butterfly until the attack was launched                            :  4  
 mean time in seconds from the beginning of a trial to when a blue tit first visited the log on the floor (first visit):  3  
 AM culculated the survival rate from raw data                                                                         :  2  
 2 cm length                                                                                                           :  1  
 (Other)                                                                                                               :  5  
   Treatment_stimulus Number_pattern   Diameter_pattern  Area_pattern   
 conspicuous: 84      Min.   : 1.000   Min.   : 1.70    Min.   :  2.27  
 eyespot    :186      1st Qu.: 2.000   1st Qu.: 4.50    1st Qu.: 15.90  
                      Median : 2.000   Median : 7.00    Median : 38.48  
                      Mean   : 2.293   Mean   :10.46    Mean   :130.34  
                      3rd Qu.: 2.000   3rd Qu.:11.48    3rd Qu.: 79.64  
                      Max.   :11.000   Max.   :31.28    Max.   :490.87  
                                                                        
 Area_background     Shape_pattern   Control_stimulus
 Min.   :  129.6   circle   : 60   bh        :  1    
 1st Qu.:  770.0   diamond  :  5   bhv       :  1    
 Median : 1040.0   eyespot  :186   bv        :  1    
 Mean   : 3040.4   lectangle: 14   bvv       :  1    
 3rd Qu.: 1444.0   others   :  4   camouflage: 10    
 Max.   :26400.0   stripe   :  1   feces     :  3    
                                   no pattern:253    
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Note_stimulus
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    : 30  
 eyespots-defensive posture vs. no eyespots-defensive posture ;models were based loosely on the late instars of Papilio canadensis and Papilio glaucus caterpillars                                                                                                                                                                                                                                                                                                                                                                                 : 12  
 single, paper models of the dingy bush brown butterfly; model which has reduced wings and eyespot diameter of approximately 6_mm. A mirror image of the butterfly was created, and a band (approx. 6_mm in length and 5_mm in width) with the same colour as the butterfly's body was added between the two images connecting the mid-section of the body. The butterfly images were then printed in colour on normal paper using a colour laser printer. The paper butterflies were cut out and dipped in paraffin wax to render them waterproof. :  9  
 solid-eyespots vs. solid-no eyespots ;models were based loosely on the late instars of Papilio canadensis and Papilio glaucus caterpillars                                                                                                                                                                                                                                                                                                                                                                                                         :  9  
 countershaded-eyespots vs. countershaded-no eyespots ;models were based loosely on the late instars of Papilio canadensis and Papilio glaucus caterpillars                                                                                                                                                                                                                                                                                                                                                                                         :  8  
 models were based loosely on the late instars of Papilio canadensis and Papilio glaucus caterpillars                                                                                                                                                                                                                                                                                                                                                                                                                                               :  7  
 (Other)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            :195  
      Type_prey                  Shape_prey        Diet_bird     Bird_sex  
 artificial:204   abstract_butterfly  :78   grainivore  :  1   both  : 50  
 real      : 66   abstract_caterpillar:68   Invertebrate:  4   female: 12  
                  abstract_prey       :58   omnivore    :122   male  : 16  
                  butterfly           :66   NA's        :143   NA's  :192  
                                                                           
                                                                           
                                                                           
     Bird_age         Tn               Cn            T_mean         
 adult   : 48   Min.   :  4.00   Min.   :  4.0   Min.   :   0.0004  
 chick   : 38   1st Qu.:  8.25   1st Qu.:  9.0   1st Qu.:   0.6021  
 juvenile:  4   Median : 20.00   Median : 20.0   Median :   7.2595  
 nestling:  1   Mean   : 42.09   Mean   : 41.9   Mean   : 100.4589  
 NA's    :179   3rd Qu.: 80.00   3rd Qu.: 80.0   3rd Qu.:  93.1250  
                Max.   :160.00   Max.   :160.0   Max.   :1180.4878  
                                                 NA's   :150        
     C_mean             MT_mean             MC_mean          
 Min.   :   0.0400   Min.   :-556.1000   Min.   :-2115.0000  
 1st Qu.:   0.3675   1st Qu.:  -1.3776   1st Qu.:   -1.1749  
 Median :   3.9950   Median :  -0.1148   Median :   -0.1552  
 Mean   :  92.8595   Mean   : -26.1327   Mean   :  -72.1179  
 3rd Qu.:  56.6866   3rd Qu.:   0.8822   3rd Qu.:    0.8949  
 Max.   :2115.0000   Max.   :   0.9996   Max.   :    0.9584  
 NA's   :150         NA's   :209         NA's   :209         
      T_sd               C_sd           T_proportion     C_proportion   
 Min.   :  0.0100   Min.   :   0.050   Min.   :0.0100   Min.   :0.0100  
 1st Qu.:  0.5278   1st Qu.:   0.229   1st Qu.:0.1715   1st Qu.:0.1046  
 Median :  4.2319   Median :   2.711   Median :0.3152   Median :0.2027  
 Mean   : 67.7837   Mean   :  64.351   Mean   :0.3531   Mean   :0.2831  
 3rd Qu.: 59.4457   3rd Qu.:  52.998   3rd Qu.:0.5000   3rd Qu.:0.4396  
 Max.   :731.8482   Max.   :1052.175   Max.   :0.9000   Max.   :0.8333  
 NA's   :147        NA's   :147        NA's   :94       NA's   :94      
      Study_design     Dataset    Mean_median  T_mean_median      
 dependent  : 65   predator:117   mea   :  1   Min.   :   0.0004  
 independent:205   prey    :153   mean  :100   1st Qu.:   0.6584  
                                  median: 19   Median :   7.2766  
                                  NA's  :150   Mean   :  89.8277  
                                               3rd Qu.:  58.0122  
                                               Max.   :1180.4878  
                                               NA's   :134        
 C_mean_median       Type_of_variance_statistic      T_se          
 Min.   :   0.0000   IQR : 19                   Min.   :  0.00000  
 1st Qu.:   0.4022   SD  :  3                   1st Qu.:  0.08208  
 Median :   4.7365   SE  : 29                   Median :  3.52000  
 Mean   :  83.4133   SE  :  1                   Mean   : 24.82370  
 3rd Qu.:  46.5061   SEM : 68                   3rd Qu.: 24.55107  
 Max.   :2115.0000   NA's:150                   Max.   :182.92683  
 NA's   :134                                    NA's   :169        
      C_se          T_variance_q1      T_variance_q3      C_variance_q1    
 Min.   :  0.0000   Min.   : 0.07075   Min.   :  0.1651   Min.   : 0.1136  
 1st Qu.:  0.0877   1st Qu.: 0.28640   1st Qu.:  0.6656   1st Qu.: 0.2594  
 Median :  1.1000   Median : 0.50720   Median :  0.9296   Median : 0.3440  
 Mean   : 25.6820   Mean   : 4.57594   Mean   : 14.0085   Mean   : 1.9774  
 3rd Qu.: 19.8747   3rd Qu.: 0.78560   3rd Qu.:  1.7799   3rd Qu.: 0.8048  
 Max.   :372.0000   Max.   :77.50000   Max.   :246.0000   Max.   :29.0000  
 NA's   :169        NA's   :251        NA's   :251        NA's   :251      
 C_variance_q3    
 Min.   : 0.2096  
 1st Qu.: 0.4340  
 Median : 0.5936  
 Mean   : 2.6702  
 3rd Qu.: 1.4151  
 Max.   :36.0000  
 NA's   :251      
                                                                                                                                          Note    
                                                                                                                                            :218  
  authors only reported total sample size (left eyespots and right eyespots), so the sample size was equally divided between the two stimuli:  8  
 Mann-Whitney U test; p < 0.002                                                                                                             :  4  
 ｱ 1 SEM                                                                                                                                    :  4  
 Calculated by AM                                                                                                                           :  3  
 (Other)                                                                                                                                    : 31  
 NA's                                                                                                                                       :  2  
      lnRR             lnRR_var             Obs_ID        Log_diameter   
 Min.   :-2.15769   Min.   :0.0008295   Min.   :  1.00   Min.   :0.5306  
 1st Qu.:-0.03772   1st Qu.:0.0111211   1st Qu.: 68.25   1st Qu.:1.5041  
 Median : 0.14784   Median :0.0234099   Median :135.50   Median :1.9459  
 Mean   : 0.26023   Mean   :0.0872373   Mean   :135.50   Mean   :2.0819  
 3rd Qu.: 0.47455   3rd Qu.:0.0732753   3rd Qu.:202.75   3rd Qu.:2.4407  
 Max.   : 4.16777   Max.   :1.3885984   Max.   :270.00   Max.   :3.4430  
                                                                         
    Log_area      Log_background  
 Min.   :0.8197   Min.   : 4.864  
 1st Qu.:2.7666   1st Qu.: 6.646  
 Median :3.6503   Median : 6.947  
 Mean   :3.8513   Mean   : 7.170  
 3rd Qu.:4.3774   3rd Qu.: 7.275  
 Max.   :6.1962   Max.   :10.181  
                                  

Final dataset for meta-analysis

I cannot attach the caption to datatable() format table. I need to use kable() format table?

Code

datatable(dt_all, caption = "Dataset for meta-analysis", 
          options = list(pageLength = 10, scrollX = TRUE))

Table 2 - Dataset for meta-analysis

Check phylogenetic relatedness

If phylogeny do not explain heterogeniety much, we will not need to consider it and the two datasets can be merged.

We used meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins. et al., 2003) and the partial I2 explained by each random factor.

Based on this, we need not to consider phylogenetic relatedness from our meta-regressions as this random factor did not explain much of the heterogeneity between effect sizes (partial I2 < 0.001%), then we can combine two datasets (predator and prey) for meta-analysis.

Code

phy_model <- function(cor_tree = vcv_tree){
            model <- rma.mv(yi = lnRR,
                            V = VCV_pred,
                            random = list(~1 | Study_ID,
                                          ~1 | Shared_control_ID,
                                          ~1 | Cohort_ID,
                                          ~1 | Obs_ID,
                                          ~1 | Bird_species,
                                          ~1 | Phylogeny),
                            R = list(Phylogeny = cor_tree),
                            test = "t",
                            method = "REML",
                            sparse = TRUE,
                            data = dat_pred)
  model
}

tree_50 <- trees[1:50]
vcv_tree_50 <- map(tree_50, ~vcv(.x, corr = TRUE))

# Run multiple meta-analyses with 50 different trees and obtain the combined result

ma_50 <- mclapply(vcv_tree_50, phy_model, mc.cores = 8) 

Code

ma_50 <- readRDS(here("Rdata", "ma_50.RDS"))

# combining the results
est_50 <- map_dbl(ma_50, ~ .x$b[[1]])
se_50 <-  map_dbl(ma_50, ~ .x$se)
df_50 <- c(rbind(est_50, se_50))

# creating an array required by pool.mi
my_array <- array(df_50, dim = c(1, 2, 50))

com_res <- round(pool.mi(my_array), 4)

Code

knitr::kable(com_res, caption = "Combined result of 50 phylogenetic trees")

Combined result of 50 phylogenetic trees

pooled.mean	pooled.total.se	se.within	se.between	relative.inc.var	frac.miss.info	CI.1	CI.2	p.value
0.1394	0.1192	0.1192	0	0	0	-0.0943	0.3731	0.2424

Code

sigma2_mod <- do.call(rbind, lapply(ma_50, function(x) x$sigma2))
sigma2_mod <- data.frame(sigma2_mod)

colnames(sigma2_mod) <- c("sigma^2.1_Study_ID", "sigma^2.2_SharedControl_ID", 
                          "sigma^2.3_Cohort_ID", "sigma^2.4_Obs_ID", 
                          "sigma^2.5_BirdSpecies", "sigma^2.6_Phylogeny")

sigma2_mod <- round(colMeans(sigma2_mod), 4)

knitr::kable(sigma2_mod, caption = "Average variance components")

Average variance components

	x
sigma^2.1_Study_ID	0.0000
sigma^2.2_SharedControl_ID	0.0923
sigma^2.3_Cohort_ID	0.1009
sigma^2.4_Obs_ID	0.5323
sigma^2.5_BirdSpecies	0.0000
sigma^2.6_Phylogeny	0.0000

We got 1000 phylogenetic trees from BirdTree.org

Only 50 trees are used as in Nakagawa & Villemereuil (2019)

2 Meta-analysis - predator dataset

We used meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins. et al., 2003) and the partial I2 explained by each random factor.

Which random effects remove in our models?

We can remove Study_ID

Code

ma_pred <- rma.mv(yi = lnRR,
                  V = VCV_pred, 
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_pred)

summary(ma_pred)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-150.0443   300.0886   310.0886   323.8565   310.6340   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0000  0.0001     18     no           Study_ID 
sigma^2.2  0.0786  0.2803     30     no  Shared_control_ID 
sigma^2.3  0.1271  0.3565     34     no          Cohort_ID 
sigma^2.4  0.4563  0.6755    117     no             Obs_ID 

Test for Heterogeneity:
Q(df = 116) = 1280.1795, p-val < .0001

Model Results:

estimate      se    tval   df    pval    ci.lb   ci.ub    
  0.1657  0.1164  1.4234  116  0.1573  -0.0649  0.3963    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_i2 <- round(i2_ml(ma_pred), 4)

pred_r2 <- round(r2_ml(ma_pred), 4)

I2_Total	I2_Study_ID	I2_Shared_control_ID	I2_Cohort_ID	I2_Obs_ID
98.8884	0	11.7364	18.9879	68.1641

R2_marginal	R2_conditional
0	0.3107

Table 3 - Heterogeneity of effect size

Code

orchard_plot(ma_pred,
              group = "Study_ID",
              xlab = "log response ratio (lnRR)", 
              angle = 45) + 
              scale_x_discrete(labels = c("Overall effect")) + 
              scale_colour_okabe_ito()+
              scale_fill_okabe_ito()

Figure 2— Estimates of overpred effect size and 95% confidence intervals

2.1 Meta-regressions: uni-moderator - predetor dataset

Pattern type

Eyespot vs. conspicuous patterns - Is there significant difference of effect size between two patterns?

Code

#normal model
pred_mr_eyespot <- rma.mv(yi = lnRR,
                    V = VCV_pred, 
                    mods = ~ Treatment_stimulus,
                    random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                    test = "t",
                    method = "REML", 
                    sparse = TRUE,
                    data = dt_pred)

summary(pred_mr_eyespot)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.9140   297.8280   307.8280   321.5527   308.3785   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0833  0.2885     30     no  Shared_control_ID 
sigma^2.2  0.1229  0.3506     34     no          Cohort_ID 
sigma^2.3  0.4600  0.6782    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1278.4557, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.6546, p-val = 0.4201

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.0809  0.3268  -0.2475  115  0.8049  -0.7283 
Treatment_stimuluseyespot    0.2812  0.3476   0.8091  115  0.4201  -0.4073 
                            ci.ub    
intrcpt                    0.5665    
Treatment_stimuluseyespot  0.9697    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_1 <- round(r2_ml(pred_mr_eyespot), 4)

R2_marginal	R2_conditional
0.0067	0.3141

Table XX. Model goodness-of-fit

Code

orchard_plot(pred_mr_eyespot,
            mod = "Treatment_stimulus",
            group = "Study_ID",
            xlab = "log response ratio (lnRR)",
            angle = 45) + 
            scale_colour_okabe_ito(order = 2:3)+
            scale_fill_okabe_ito(order = 2:3)

Figure 3— Eyespot vs. conspicuous patterns

Code

# intercept-removed model
pred_mr_eyespot1 <- rma.mv(yi = lnRR,
                      V = VCV_pred, 
                      mods = ~ Treatment_stimulus -1,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                                test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_pred)

summary(pred_mr_eyespot1)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.9140   297.8280   307.8280   321.5527   308.3785   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0833  0.2885     30     no  Shared_control_ID 
sigma^2.2  0.1229  0.3506     34     no          Cohort_ID 
sigma^2.3  0.4600  0.6782    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1278.4557, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 115) = 1.3386, p-val = 0.2663

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
Treatment_stimulusconspicuous   -0.0809  0.3268  -0.2475  115  0.8049  -0.7283 
Treatment_stimuluseyespot        0.2003  0.1242   1.6127  115  0.1096  -0.0457 
                                ci.ub    
Treatment_stimulusconspicuous  0.5665    
Treatment_stimuluseyespot      0.4463    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pattern diameter (mm)

Code

pred_mr_diameter <- rma.mv(yi = lnRR,
                      V = VCV_pred, 
                      mods = ~ Log_diameter,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_pred)

summary(pred_mr_diameter)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.6526   297.3051   307.3051   321.0298   307.8556   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0829  0.2879     30     no  Shared_control_ID 
sigma^2.2  0.1316  0.3628     34     no          Cohort_ID 
sigma^2.3  0.4566  0.6758    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1269.1592, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.3781, p-val = 0.5398

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt        -0.0240  0.3296  -0.0729  115  0.9420  -0.6768  0.6288    
Log_diameter    0.0880  0.1432   0.6149  115  0.5398  -0.1955  0.3716    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_2 <- round(r2_ml(pred_mr_diameter), 4)

R2_marginal	R2_conditional
0.0085	0.3254

Table XX. Model goodness-of-fit

Code

bubble_plot(pred_mr_diameter,
            mod = "Log_diameter",
            group = "Study_ID",
            xlab = "Log-transformed diameter")

Pattern area (mm²)

Code

pred_mr_area <- rma.mv(yi = lnRR,
                  V = VCV_pred, 
                  mods = ~ Log_area,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_pred)

summary(pred_mr_area)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.5635   297.1269   307.1269   320.8516   307.6774   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0816  0.2857     30     no  Shared_control_ID 
sigma^2.2  0.1298  0.3603     34     no          Cohort_ID 
sigma^2.3  0.4571  0.6761    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1269.5060, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.5212, p-val = 0.4718

Model Results:

          estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt    -0.0476  0.3175  -0.1498  115  0.8812  -0.6766  0.5814    
Log_area    0.0529  0.0732   0.7219  115  0.4718  -0.0922  0.1979    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_3 <- round(r2_ml(pred_mr_area), 4)

R2_marginal	R2_conditional
0.0122	0.3245

Table XX. Model goodness-of-fit

Code

bubble_plot(pred_mr_area,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Pattern number

Code

pred_mr_num <- rma.mv(yi = lnRR,
                  V = VCV_pred,
                  mods = ~ Number_pattern,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 |Obs_ID),                      
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_pred)

summary(pred_mr_num)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-149.1039   298.2077   308.2077   321.9324   308.7582   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0838  0.2894     30     no  Shared_control_ID 
sigma^2.2  0.1289  0.3590     34     no          Cohort_ID 
sigma^2.3  0.4602  0.6784    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1277.6943, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.0197, p-val = 0.8885

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt           0.2002  0.2728   0.7338  115  0.4645  -0.3402  0.7405    
Number_pattern   -0.0161  0.1145  -0.1405  115  0.8885  -0.2430  0.2108    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_4 <- round(r2_ml(pred_mr_num), 4)

R2_marginal	R2_conditional
3e-04	0.3162

Table XX. Model goodness-of-fit

Code

bubble_plot(pred_mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") 

Prey material type

Our dataset includes two material types of prey: live and artificial. Is there significant difference of effect size between two stimuli?

Code

# normal model
pred_mr_prey_type <- rma.mv(yi = lnRR,
                  V = VCV_pred, 
                  mods = ~ Type_prey,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID), 
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_pred)

summary(pred_mr_prey_type)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.7653   297.5306   307.5306   321.2553   308.0810   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0793  0.2816     30     no  Shared_control_ID 
sigma^2.2  0.1530  0.3912     34     no          Cohort_ID 
sigma^2.3  0.4487  0.6699    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1232.0173, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.4715, p-val = 0.4937

Model Results:

               estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt          0.0994  0.1519  0.6543  115  0.5142  -0.2015  0.4002    
Type_preyreal    0.1699  0.2474  0.6866  115  0.4937  -0.3202  0.6600    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_5 <- round(r2_ml(pred_mr_prey_type), 4)

R2_marginal	R2_conditional
0.0079	0.3463

Table XX. Model goodness-of-fit

Code

orchard_plot(pred_mr_prey_type,
              mod = "Type_prey",
              group = "Study_ID",
              xlab = "Type of prey",
              angle = 45) +
            scale_colour_okabe_ito(order = 4:5)+
            scale_fill_okabe_ito(order = 4:5)

Code

# intercept-removed model
pred_mr_prey_type1 <- rma.mv(yi = lnRR,
                            V = VCV_pred, 
                            mods = ~ Type_prey -1,
                            random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID), 
                            test = "t",
                            method = "REML", 
                            sparse = TRUE,
                            data = dt_pred)

summary(pred_mr_prey_type1)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.7653   297.5306   307.5306   321.2553   308.0810   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0793  0.2816     30     no  Shared_control_ID 
sigma^2.2  0.1530  0.3912     34     no          Cohort_ID 
sigma^2.3  0.4487  0.6699    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1232.0173, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 115) = 1.1641, p-val = 0.3159

Model Results:

                     estimate      se    tval   df    pval    ci.lb   ci.ub    
Type_preyartificial    0.0994  0.1519  0.6543  115  0.5142  -0.2015  0.4002    
Type_preyreal          0.2693  0.1953  1.3785  115  0.1707  -0.1177  0.6562    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prey shape type

Our dataset includes 4 types of prey type: real butterfly, abstract butterfly, abstract caterpillar, and abstract prey. Is there significant difference of effect size between two stimuli?

Code

pred_mr_prey_shape <- rma.mv(yi = lnRR,
                        V = VCV_pred, 
                        mods = ~ Shape_prey -1,
                        random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),,
                                      test = "t",
                                      method = "REML", 
                                      sparse = TRUE,
                                      data = dt_pred)

summary(pred_mr_prey_shape)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-146.3709   292.7417   306.7417   325.8334   307.8084   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0998  0.3159     30     no  Shared_control_ID 
sigma^2.2  0.1719  0.4146     34     no          Cohort_ID 
sigma^2.3  0.4475  0.6689    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 113) = 1227.2585, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 113) = 0.5653, p-val = 0.6883

Model Results:

                                estimate      se    tval   df    pval    ci.lb 
Shape_preyabstract_butterfly      0.1090  0.3403  0.3203  113  0.7493  -0.5652 
Shape_preyabstract_caterpillar    0.0181  0.2757  0.0656  113  0.9478  -0.5280 
Shape_preyabstract_prey           0.1425  0.2407  0.5920  113  0.5551  -0.3344 
Shape_preybutterfly               0.2724  0.2028  1.3431  113  0.1819  -0.1294 
                                 ci.ub    
Shape_preyabstract_butterfly    0.7832    
Shape_preyabstract_caterpillar  0.5642    
Shape_preyabstract_prey         0.6194    
Shape_preybutterfly             0.6742    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

pred_r2_6 <- round(r2_ml(pred_mr_prey_shape), 4)

R2_marginal	R2_conditional
0.0088	0.3832

Table XX. Model goodness-of-fit

Code

orchard_plot(pred_mr_prey_shape,
              mod = "Shape_prey",
              group = "Study_ID",
              xlab = "Shape of prey",
              angle = 45) +
              scale_colour_okabe_ito(order = 6:9)+
              scale_fill_okabe_ito(order = 6:9)

Figure 4— Effect size and prey shape

2.2 Correlation visualisation and choose moderators for multi-moderator meta-regression in predator dataset

Before we run multi-moderator meta-regressions, we need to consider the correlation between moderators. Area, diameter and background seem to be correlated. We visualised the correlation between these variables.

Visualisation of correlation

Code

# correlation between continuous variables
pred_corr_cont <- round(cor(dt_pred[, c("Diameter_pattern", "Area_pattern", "Number_pattern", "Area_background")]),2)

pred_p1 <- ggcorrplot(pred_corr_cont, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(a) Continuous variables")

pred_corr_cont_log <- round(cor(dt_pred[, c("Log_diameter", "Log_area", "Number_pattern", "Log_background")]),2)

pred_p2 <- ggcorrplot(pred_corr_cont_log, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(b) Log-transormed continuous variables")

pred_p1 + pred_p2 +plot_layout(guides = 'collect')

Figure 5— Correlation between coninuous moderators

Code

summary(dt_pred$Treatment_stimulus)

conspicuous     eyespot 
          7         110 

Code

pred_dat1 <- dt_pred %>%
  dplyr::select("Treatment_stimulus", "Type_prey", "Shape_prey")

pred_corr_cat <- GKtauDataframe(pred_dat1)
plot(pred_corr_cat)

Figure 6— Correlation between categorical variables

We should not include “Shape_prey” and “Type_prey” in the model at the same time. Because they are correlated.

Choose moderators

We used model R2 values to find better model and modelator VIF values to check multicollinearity between moderators. Higher R2 indicates better model and VIF > 2 indicates multicollinearity.

2.2.0.1 R2

Code

pred_r2_area <- rma.mv(yi = lnRR,
                  V = VCV_pred,
                  mods = ~ Log_area + Number_pattern,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML",
                  sparse = TRUE,
                  data = dt_pred)

pred_r2_diameter <- rma.mv(yi = lnRR,
                      V = VCV_pred,
                      mods = ~ Log_diameter + Number_pattern,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_pred)

# check r2 values
r2_ml(pred_r2_area)

   R2_marginal R2_conditional 
    0.01211581     0.33044840 

Code

r2_ml(pred_r2_diameter)

   R2_marginal R2_conditional 
   0.008825168    0.330370102 

It seems area is a better predictor than diameter .

2.3 Meta-regressions: multi-moderator

All potential moderators (full model)

Include all moderators in the model.

Code

pred_mr_all <- rma.mv(yi = lnRR,
                V = VCV_pred, 
                mods = ~ Treatment_stimulus + Log_area 
                + Number_pattern + Type_prey,
                random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_pred)

summary(pred_mr_all)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-144.9280   289.8560   305.8560   327.6040   307.2541   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0859  0.2931     30     no  Shared_control_ID 
sigma^2.2  0.1358  0.3685     34     no          Cohort_ID 
sigma^2.3  0.4645  0.6816    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 112) = 1225.0519, p-val < .0001

Test of Moderators (coefficients 2:5):
F(df1 = 4, df2 = 112) = 0.5654, p-val = 0.6883

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.3783  0.5230  -0.7234  112  0.4709  -1.4146 
Treatment_stimuluseyespot    0.2844  0.3530   0.8056  112  0.4222  -0.4151 
Log_area                     0.0815  0.0789   1.0329  112  0.3039  -0.0748 
Number_pattern              -0.0657  0.1246  -0.5271  112  0.5992  -0.3127 
Type_preyreal                0.2876  0.2785   1.0328  112  0.3039  -0.2642 
                            ci.ub    
intrcpt                    0.6579    
Treatment_stimuluseyespot  0.9839    
Log_area                   0.2378    
Number_pattern             0.1813    
Type_preyreal              0.8395    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(pred_mr_all)

   R2_marginal R2_conditional 
    0.03194282     0.34468854 

Moderators associated with conspicuousness

Only include moderators related to conspicuousness (pattern area and number of patterns) in the model.

Code

pred_mr_cons <- rma.mv(yi = lnRR,
                V = VCV_pred, 
                mods = ~ Log_area + Number_pattern,
                random = list(~1 | Shared_control_ID,
                              ~1 | Cohort_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_pred)

summary(pred_mr_cons)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-147.6172   295.2344   307.2344   323.6516   308.0195   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0865  0.2942     30     no  Shared_control_ID 
sigma^2.2  0.1326  0.3642     34     no          Cohort_ID 
sigma^2.3  0.4609  0.6789    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 114) = 1268.2068, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 114) = 0.2623, p-val = 0.7698

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          -0.0122  0.4025  -0.0304  114  0.9758  -0.8096  0.7852    
Log_area          0.0529  0.0742   0.7126  114  0.4775  -0.0941  0.1998    
Number_pattern   -0.0164  0.1152  -0.1425  114  0.8869  -0.2447  0.2118    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(pred_mr_cons)

   R2_marginal R2_conditional 
    0.01211581     0.33044840 

Code

bubble_plot(pred_mr_cons,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Code

pred_mr_cons_1 <- rma.mv(yi = lnRR,
                V = VCV_pred, 
                mods = ~ Log_area + Number_pattern + Treatment_stimulus,
                random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_pred)

summary(pred_mr_cons_1)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-146.4650   292.9300   306.9300   326.0217   307.9967   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0903  0.3005     30     no  Shared_control_ID 
sigma^2.2  0.1320  0.3633     34     no          Cohort_ID 
sigma^2.3  0.4635  0.6808    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 113) = 1264.4653, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 113) = 0.3975, p-val = 0.7550

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.2715  0.5129  -0.5293  113  0.5976  -1.2877 
Log_area                     0.0548  0.0746   0.7340  113  0.4644  -0.0931 
Number_pattern              -0.0177  0.1157  -0.1530  113  0.8787  -0.2470 
Treatment_stimuluseyespot    0.2902  0.3529   0.8222  113  0.4127  -0.4090 
                            ci.ub    
intrcpt                    0.7447    
Log_area                   0.2027    
Number_pattern             0.2116    
Treatment_stimuluseyespot  0.9894    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(pred_mr_cons_1)

   R2_marginal R2_conditional 
    0.02106248     0.33834441 

Moderators associated with prey in predator dataset

Code

pred_mr_pre <- rma.mv(yi = lnRR,
                V = VCV_pred, 
                mods = ~ Treatment_stimulus + Type_prey,
                random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_pred)

summary(pred_mr_pre)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-147.6573   295.3146   307.3146   323.7318   308.0997   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0844  0.2905     30     no  Shared_control_ID 
sigma^2.2  0.1473  0.3838     34     no          Cohort_ID 
sigma^2.3  0.4534  0.6733    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 114) = 1230.6103, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 114) = 0.5265, p-val = 0.5921

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.1347  0.3416  -0.3943  114  0.6941  -0.8115 
Treatment_stimuluseyespot    0.2714  0.3530   0.7689  114  0.4435  -0.4278 
Type_preyreal                0.1622  0.2481   0.6536  114  0.5147  -0.3294 
                            ci.ub    
intrcpt                    0.5421    
Treatment_stimuluseyespot  0.9706    
Type_preyreal              0.6538    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(pred_mr_pre)

   R2_marginal R2_conditional 
    0.01177801     0.34603083 

2.4 Publication bias

Funnel plot

Code

# funnel plot - standard error
funnel(ma_pred, yaxis = "sei", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)" )

# funnel plot - inverse of standard error
pred_f2 <- funnel(ma_pred, yaxis = "seinv", 
      xlab = "Standarized residuals",
      ylab = "Precision (inverse of SE)",  col = c(alpha("orange", 0.5)))

Figure 7— Effect size and its standard error

Figure 8— Effect size and its inverse standard error

Figure 9— Funnel plot of effect size and its standard error

Egger’s test

Code

pred_df_bias <- dt_pred %>% mutate(sqrt_inv_e_n = sqrt((Cn + Tn)/(Cn * Tn)))

pred_bias_model <- rma.mv(yi = lnRR,
                      V = lnRR_var, 
                      mods = ~ 1 + sqrt_inv_e_n,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = pred_df_bias)

summary(pred_bias_model)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-149.0844   298.1687   308.1687   321.8934   308.7192   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0804  0.2835     30     no  Shared_control_ID 
sigma^2.2  0.1285  0.3585     34     no          Cohort_ID 
sigma^2.3  0.4601  0.6783    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1249.3192, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.0592, p-val = 0.8082

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt         0.2525  0.3752   0.6729  115  0.5023  -0.4907  0.9957    
sqrt_inv_e_n   -0.2161  0.8884  -0.2433  115  0.8082  -1.9758  1.5435    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(pred_bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Figure 10— Egger’s test of publication bias

Decline effect

Code

pred_year_model <- rma.mv(yi = lnRR,
                      V = VCV_pred, 
                      mods = ~  1 + Year,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = pred_df_bias)

summary(pred_year_model)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.7939   297.5878   307.5878   321.3125   308.1383   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0843  0.2904     30     no  Shared_control_ID 
sigma^2.2  0.1359  0.3687     34     no          Cohort_ID 
sigma^2.3  0.4563  0.6755    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1254.3394, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.0430, p-val = 0.8361

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt    4.3742  20.3031   0.2154  115  0.8298  -35.8423  44.5907    
Year      -0.0021   0.0101  -0.2073  115  0.8361   -0.0221   0.0179    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(pred_year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

Figure 11— Decline effect

multi-moderator model

Code

pred_bias_all <- rma.mv(yi = lnRR,
                      V = VCV_pred, 
                      mods = ~ 1 + sqrt_inv_e_n + Year,
                      random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = pred_df_bias)

summary(pred_bias_all)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-147.7933   295.5867   307.5867   324.0038   308.3717   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0852  0.2919     30     no  Shared_control_ID 
sigma^2.2  0.1372  0.3704     34     no          Cohort_ID 
sigma^2.3  0.4606  0.6786    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 114) = 1225.6345, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 114) = 0.0569, p-val = 0.9447

Model Results:

              estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt         5.6563  20.9515   0.2700  114  0.7877  -35.8483  47.1610    
sqrt_inv_e_n   -0.2454   0.9196  -0.2668  114  0.7901   -2.0670   1.5763    
Year           -0.0027   0.0104  -0.2584  114  0.7966   -0.0233   0.0179    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.5 Additional analyses

2.5.1 Background area (mm²)

Code

pred_mr_background <- rma.mv(yi = lnRR,
                        V = VCV_pred, 
                        mods = ~ Log_background,
                  random = list(~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                        test = "t",
                        method = "REML",
                        sparse = TRUE,
                        data = dt_pred)

summary(pred_mr_background)

Multivariate Meta-Analysis Model (k = 117; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-148.6517   297.3033   307.3033   321.0280   307.8538   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0960  0.3098     30     no  Shared_control_ID 
sigma^2.2  0.1420  0.3768     34     no          Cohort_ID 
sigma^2.3  0.4490  0.6701    117     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 115) = 1270.4543, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 115) = 0.2344, p-val = 0.6292

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt           0.4975  0.6997   0.7110  115  0.4785  -0.8885  1.8834    
Log_background   -0.0446  0.0922  -0.4842  115  0.6292  -0.2272  0.1380    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(pred_mr_background,
            mod = "Log_background",
            group = "Study_ID",
            xlab = "Log-transformed backgroud")

Figure 12— Effect size and background area

2.6 Figure galary

Discrete moderators

Code

# overall effect
pred_p1 <- orchard_plot(ma_pred,
                    group = "Study_ID",
                    xlab = "log response ratio (lnRR)", 
                    angle = 45,
                    legend.pos = "bottom.right") + 
                    scale_x_discrete(labels = c("Overall effect")) + 
                    theme(legend.direction = "horizontal",
                    legend.text = element_text(size = 6)) +
                    
                    scale_colour_okabe_ito()+
                    scale_fill_okabe_ito()

#eyespot vs. conspicuous patterns
pred_p2 <- orchard_plot(pred_mr_eyespot,
                  mod = "Treatment_stimulus",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6)) +
                  scale_colour_okabe_ito(order = 2:3)+
                  scale_fill_okabe_ito(order = 2:3)
                  
# prey type - real vs. artificial
pred_p3 <- orchard_plot(pred_mr_prey_type,
                  mod = "Type_prey",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6))+ 
                  scale_colour_okabe_ito(order = 4:5) +
                  scale_fill_okabe_ito(order = 4:5)

# combine
pred_patch1 <- pred_p1 | (pred_p2 / pred_p3)
pred_patch1 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig1.pdf", width = 10, height = 5, dpi = 450)

Figure 13— Discrete moderators

Continuous moderators

Code

# these figs were created by multi meta-regression model results
# log-transformed area
pred_p4 <- bubble_plot(pred_mr_area,
                  mod = "Log_area",
                  group = "Study_ID",
                  xlab = "Log-transformed area") +
                  scale_x_continuous(breaks = seq(0,7,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# number of patterns
pred_p5 <- bubble_plot(pred_mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") +
                  scale_x_continuous(breaks = seq(0,11,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# combine
pred_p4 / pred_p5 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig2_multi.pdf", width = 10, height = 10, dpi = 450)

Figure 14— Continuous moderators

Publication bias

Code

# funnel plot

# pdf("fig3.pdf")

funnel(ma_pred, yaxis = "seinv", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)",
      xlim = c(-4.0, 4.5), 
      ylim = c(0.01, 80.0), 
      col = c(alpha("mediumvioletred", 0.5)),
      cex = 0.7)

# dev.off()

# Egger's test and decline effect
pred_p7 <- bubble_plot(pred_bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size") +
            scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5)) 

pred_p8 <- bubble_plot(pred_year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication") +  scale_y_continuous(breaks = seq(-2.5, 4.0, 1.5)) 

# combine plots created by bubble_plot()
pred_pub <- pred_p7 / pred_p8
pred_pub +  plot_annotation(tag_levels = 'a') 

# output figure as pdf
# ggsave("fig4.pdf", width = 10, height = 10, dpi = 450)

Figure 15— Funnel plot

Figure 16— Egger’s test and Decline effect

Figure 17— Publication bias

3 Meta-analysis - prey dataset

We used meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins. et al., 2003) and the partial I2 explained by each random factor.

Code

ma_prey <- rma.mv(yi = lnRR,
                  V = VCV_prey, 
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_prey)

summary(ma_prey)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-51.2439  102.4878  112.4878  127.6072  112.8988   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0787  0.2806     16     no           Study_ID 
sigma^2.2  0.0246  0.1568     59     no  Shared_control_ID 
sigma^2.3  0.0027  0.0522    130     no          Cohort_ID 
sigma^2.4  0.0298  0.1727    153     no             Obs_ID 

Test for Heterogeneity:
Q(df = 152) = 905.4586, p-val < .0001

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub     
  0.2156  0.0818  2.6345  152  0.0093  0.0539  0.3773  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_i2 <- round(i2_ml(ma_prey), 4)

prey_r2 <- round(r2_ml(ma_prey), 4)

I2_Total	I2_Study_ID	I2_Shared_control_ID	I2_Cohort_ID	I2_Obs_ID
90.2479	52.2906	16.3329	1.8132	19.8112

R2_marginal	R2_conditional
0	0.7805

Table 3 - Heterogeneity of effect size

Code

orchard_plot(ma_prey,
              group = "Study_ID",
              xlab = "log response ratio (lnRR)", 
              angle = 45) + 
              scale_x_discrete(labels = c("Overall effect")) + 
              scale_colour_okabe_ito()+
              scale_fill_okabe_ito()

Figure 18— Estimates of overprey effect size and 95% confidence intervals

3.1 Meta-regressions: uni-moderator - preyetor dataset

Pattern type

Eyespot vs. conspicuous patterns - Is there significant difference of effect size between two patterns?

Code

#normal model
prey_mr_eyespot <- rma.mv(yi = lnRR,
                    V = VCV_prey, 
                    mods = ~ Treatment_stimulus,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                    test = "t",
                    method = "REML", 
                    sparse = TRUE,
                    data = dt_prey)

summary(prey_mr_eyespot)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-50.1903  100.3806  112.3806  130.4843  112.9640   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0848  0.2912     16     no           Study_ID 
sigma^2.2  0.0285  0.1687     59     no  Shared_control_ID 
sigma^2.3  0.0008  0.0291    130     no          Cohort_ID 
sigma^2.4  0.0290  0.1703    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 896.1286, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 151) = 1.9624, p-val = 0.1633

Model Results:

                           estimate      se    tval   df    pval    ci.lb 
intrcpt                      0.1582  0.0943  1.6773  151  0.0956  -0.0282 
Treatment_stimuluseyespot    0.1036  0.0740  1.4009  151  0.1633  -0.0425 
                            ci.ub    
intrcpt                    0.3446  . 
Treatment_stimuluseyespot  0.2498    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_1 <- round(r2_ml(prey_mr_eyespot), 4)

R2_marginal	R2_conditional
0.0185	0.8011

Table XX. Model goodness-of-fit

Code

orchard_plot(prey_mr_eyespot,
            mod = "Treatment_stimulus",
            group = "Study_ID",
            xlab = "log response ratio (lnRR)",
            angle = 45) + 
            scale_colour_okabe_ito(order = 2:3)+
            scale_fill_okabe_ito(order = 2:3)

Figure 19— Eyespot vs. conspicuous patterns

Code

# intercept-removed model
prey_mr_eyespot1 <- rma.mv(yi = lnRR,
                      V = VCV_prey, 
                      mods = ~ Treatment_stimulus -1,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_prey)

summary(prey_mr_eyespot1)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-50.1903  100.3806  112.3806  130.4843  112.9640   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0848  0.2912     16     no           Study_ID 
sigma^2.2  0.0285  0.1687     59     no  Shared_control_ID 
sigma^2.3  0.0008  0.0291    130     no          Cohort_ID 
sigma^2.4  0.0290  0.1703    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 896.1286, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 151) = 4.2465, p-val = 0.0161

Model Results:

                               estimate      se    tval   df    pval    ci.lb 
Treatment_stimulusconspicuous    0.1582  0.0943  1.6773  151  0.0956  -0.0282 
Treatment_stimuluseyespot        0.2618  0.0907  2.8875  151  0.0045   0.0827 
                                ci.ub     
Treatment_stimulusconspicuous  0.3446   . 
Treatment_stimuluseyespot      0.4410  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pattern diameter (mm)

Code

prey_mr_diameter <- rma.mv(yi = lnRR,
                      V = VCV_prey, 
                      mods = ~ Log_diameter,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_prey)

summary(prey_mr_diameter)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-37.3428   74.6856   86.6856  104.7893   87.2689   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0306  0.1750     16     no           Study_ID 
sigma^2.2  0.0169  0.1299     59     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    130     no          Cohort_ID 
sigma^2.4  0.0283  0.1682    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 646.3148, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 151) = 32.7213, p-val < .0001

Model Results:

              estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt        -0.5415  0.1415  -3.8280  151  0.0002  -0.8210  -0.2620  *** 
Log_diameter    0.4228  0.0739   5.7203  151  <.0001   0.2768   0.5689  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_2 <- round(r2_ml(prey_mr_diameter), 4)

R2_marginal	R2_conditional
0.3144	0.744

Table XX. Model goodness-of-fit

Code

bubble_plot(prey_mr_diameter,
            mod = "Log_diameter",
            group = "Study_ID",
            xlab = "Log-transformed diameter")

Pattern area (mm²)

Code

prey_mr_area <- rma.mv(yi = lnRR,
                  V = VCV_prey, 
                  mods = ~ Log_area,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_prey)

summary(prey_mr_area)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-32.7267   65.4534   77.4534   95.5571   78.0367   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0310  0.1761     16     no           Study_ID 
sigma^2.2  0.0158  0.1258     59     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    130     no          Cohort_ID 
sigma^2.4  0.0253  0.1590    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 638.1055, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 151) = 43.3344, p-val < .0001

Model Results:

          estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt    -0.6071  0.1348  -4.5027  151  <.0001  -0.8735  -0.3407  *** 
Log_area    0.2562  0.0389   6.5829  151  <.0001   0.1793   0.3330  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_3 <- round(r2_ml(prey_mr_area), 4)

R2_marginal	R2_conditional
0.3779	0.7819

Table XX. Model goodness-of-fit

Code

bubble_plot(prey_mr_area,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Pattern number

Code

prey_mr_num <- rma.mv(yi = lnRR,
                  V = VCV_prey,
                  mods = ~ Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_prey)

summary(prey_mr_num)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-39.9363   79.8726   91.8726  109.9762   92.4559   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0820  0.2864     16     no           Study_ID 
sigma^2.2  0.0116  0.1079     59     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    130     no          Cohort_ID 
sigma^2.4  0.0267  0.1633    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 827.8470, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 151) = 25.9753, p-val < .0001

Model Results:

                estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt           0.4748  0.0952   4.9890  151  <.0001   0.2867   0.6628  *** 
Number_pattern   -0.0794  0.0156  -5.0966  151  <.0001  -0.1102  -0.0486  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_4 <- round(r2_ml(prey_mr_num), 4)

R2_marginal	R2_conditional
0.1321	0.8077

Table XX. Model goodness-of-fit

Code

bubble_plot(prey_mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") 

Prey material type

Our dataset includes two material types of prey: live and artificial. Is there significant difference of effect size between two stimuli?

Code

# normal model
prey_mr_prey_type <- rma.mv(yi = lnRR,
                  V = VCV_prey, 
                  mods = ~ Type_prey,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_prey)

summary(prey_mr_prey_type)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-50.0777  100.1555  112.1555  130.2592  112.7388   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0817  0.2858     16     no           Study_ID 
sigma^2.2  0.0249  0.1578     59     no  Shared_control_ID 
sigma^2.3  0.0022  0.0468    130     no          Cohort_ID 
sigma^2.4  0.0301  0.1736    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 898.0132, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 151) = 0.6965, p-val = 0.4053

Model Results:

               estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          0.2702  0.1054   2.5640  151  0.0113   0.0620  0.4784  * 
Type_preyreal   -0.1429  0.1712  -0.8346  151  0.4053  -0.4812  0.1954    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_5 <- round(r2_ml(prey_mr_prey_type), 4)

R2_marginal	R2_conditional
0.0264	0.7887

Table XX. Model goodness-of-fit

Code

orchard_plot(prey_mr_prey_type,
              mod = "Type_prey",
              group = "Study_ID",
              xlab = "Type of prey",
              angle = 45) +
            scale_colour_okabe_ito(order = 4:5)+
            scale_fill_okabe_ito(order = 4:5)

Code

# intercept-removed model
prey_mr_prey_type1 <- rma.mv(yi = lnRR,
                            V = VCV_prey, 
                            mods = ~ Type_prey -1,
                            random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                            test = "t",
                            method = "REML", 
                            sparse = TRUE,
                            data = dt_prey)

summary(prey_mr_prey_type1)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-50.0777  100.1555  112.1555  130.2592  112.7388   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0817  0.2858     16     no           Study_ID 
sigma^2.2  0.0249  0.1578     59     no  Shared_control_ID 
sigma^2.3  0.0022  0.0468    130     no          Cohort_ID 
sigma^2.4  0.0301  0.1736    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 898.0132, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 151) = 3.7322, p-val = 0.0262

Model Results:

                     estimate      se    tval   df    pval    ci.lb   ci.ub    
Type_preyartificial    0.2702  0.1054  2.5640  151  0.0113   0.0620  0.4784  * 
Type_preyreal          0.1273  0.1349  0.9436  151  0.3469  -0.1393  0.3939    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prey shape type

Our dataset includes 4 types of prey type: real butterfly, abstract butterfly, abstract caterpillar, and abstract prey. Is there significant difference of effect size between two stimuli?

Code

prey_mr_prey_shape <- rma.mv(yi = lnRR,
                        V = VCV_prey, 
                        mods = ~ Shape_prey -1,
                        random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                        test = "t",
                        method = "REML", 
                        sparse = TRUE,
                        data = dt_prey)

summary(prey_mr_prey_shape)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-47.3536   94.7071  108.7071  129.7816  109.4959   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0603  0.2455     16     no           Study_ID 
sigma^2.2  0.0267  0.1634     59     no  Shared_control_ID 
sigma^2.3  0.0016  0.0399    130     no          Cohort_ID 
sigma^2.4  0.0301  0.1736    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 150) = 769.1828, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 150) = 4.2577, p-val = 0.0064

Model Results:

                                estimate      se    tval   df    pval    ci.lb 
Shape_preyabstract_butterfly      0.4034  0.1176  3.4299  150  0.0008   0.1710 
Shape_preyabstract_caterpillar    0.0258  0.1536  0.1682  150  0.8666  -0.2777 
Shape_preybutterfly               0.1215  0.1227  0.9901  150  0.3237  -0.1209 
                                 ci.ub      
Shape_preyabstract_butterfly    0.6358  *** 
Shape_preyabstract_caterpillar  0.3294      
Shape_preybutterfly             0.3638      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

prey_r2_6 <- round(r2_ml(prey_mr_prey_shape), 4)

R2_marginal	R2_conditional
0.1977	0.7962

Table XX. Model goodness-of-fit

Code

orchard_plot(prey_mr_prey_shape,
              mod = "Shape_prey",
              group = "Study_ID",
              xlab = "Shape of prey",
              angle = 45) +
              scale_colour_okabe_ito(order = 6:9)+
              scale_fill_okabe_ito(order = 6:9)

Figure 20— Effect size and prey shape

3.2 Correlation visualisation and choose moderators for multi-moderator meta-regression in preyator dataset

Before we run multi-moderator meta-regressions, we need to consider the correlation between moderators. Area, diameter and background seem to be correlated. We visualised the correlation between these variables.

Visualisation of correlation

Code

# correlation between continuous variables
prey_corr_cont <- round(cor(dt_prey[, c("Diameter_pattern", "Area_pattern", "Number_pattern", "Area_background")]),2)

prey_p1 <- ggcorrplot(prey_corr_cont, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(a) Continuous variables")

prey_corr_cont_log <- round(cor(dt_prey[, c("Log_diameter", "Log_area", "Number_pattern", "Log_background")]),2)

prey_p2 <- ggcorrplot(prey_corr_cont_log, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(b) Log-transormed continuous variables")

prey_p1 + prey_p2 +plot_layout(guides = 'collect')

Figure 21— Correlation between coninuous moderators

Code

prey_dat1 <- dt_prey %>%
  dplyr::select("Treatment_stimulus", "Type_prey", "Shape_prey")

prey_corr_cat <- GKtauDataframe(prey_dat1)
plot(prey_corr_cat)

Figure 22— Correlation between categorical variables

We should not include “Shape_prey” and “Type_prey” in the model at the same time. Because they are correlated.

Choose moderators

We used model R2 values to find better model and modelator VIF values to check multicollinearity between moderators. Higher R2 indicates better model and VIF > 2 indicates multicollinearity.

3.2.0.1 R2

Code

prey_r2_area <- rma.mv(yi = lnRR,
                  V = VCV_prey,
                  mods = ~ Log_area + Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML",
                  sparse = TRUE,
                  data = dt_prey)

prey_r2_diameter <- rma.mv(yi = lnRR,
                      V = VCV_prey,
                      mods = ~ Log_diameter + Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_prey)

# check r2 values
r2_ml(prey_r2_area)

   R2_marginal R2_conditional 
     0.3411669      0.7941831 

Code

r2_ml(prey_r2_diameter)

   R2_marginal R2_conditional 
     0.2751123      0.7751117 

It seems area is a better preyictor than diameter .

3.3 Meta-regressions: multi-moderator

All potential moderators (full model)

Include all moderators in the model.

Code

prey_mr_all <- rma.mv(yi = lnRR,
                V = VCV_prey, 
                mods = ~ Treatment_stimulus + Log_area 
                + Number_pattern + Type_prey,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_prey)

summary(prey_mr_all)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-29.3878   58.7756   76.7756  103.7505   78.0799   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0481  0.2194     16     no           Study_ID 
sigma^2.2  0.0139  0.1177     59     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    130     no          Cohort_ID 
sigma^2.4  0.0224  0.1497    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 148) = 604.0858, p-val < .0001

Test of Moderators (coefficients 2:5):
F(df1 = 4, df2 = 148) = 12.2472, p-val < .0001

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.4012  0.2036  -1.9704  148  0.0507  -0.8035 
Treatment_stimuluseyespot    0.1058  0.0663   1.5944  148  0.1130  -0.0253 
Log_area                     0.2069  0.0478   4.3263  148  <.0001   0.1124 
Number_pattern              -0.0384  0.0187  -2.0565  148  0.0415  -0.0754 
Type_preyreal                0.0619  0.1469   0.4212  148  0.6742  -0.2284 
                             ci.ub      
intrcpt                     0.0012    . 
Treatment_stimuluseyespot   0.2369      
Log_area                    0.3014  *** 
Number_pattern             -0.0015    * 
Type_preyreal               0.3521      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(prey_mr_all)

   R2_marginal R2_conditional 
     0.3168985      0.8186497 

Moderators associated with conspicuousness

Only include moderators related to conspicuousness (pattern area and number of patterns) in the model.

Code

prey_mr_cons <- rma.mv(yi = lnRR,
                V = VCV_prey, 
                mods = ~ Log_area + Number_pattern,
                random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_prey)

summary(prey_mr_cons)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-31.2407   62.4815   76.4815   97.5559   77.2702   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0413  0.2033     16     no           Study_ID 
sigma^2.2  0.0124  0.1115     59     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    130     no          Cohort_ID 
sigma^2.4  0.0244  0.1563    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 150) = 627.6328, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 150) = 23.2412, p-val < .0001

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub      
intrcpt          -0.3367  0.1945  -1.7310  150  0.0855  -0.7210  0.0476    . 
Log_area          0.2063  0.0470   4.3904  150  <.0001   0.1134  0.2991  *** 
Number_pattern   -0.0340  0.0177  -1.9193  150  0.0568  -0.0690  0.0010    . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(prey_mr_cons)

   R2_marginal R2_conditional 
     0.3411669      0.7941831 

Code

bubble_plot(prey_mr_cons,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Code

prey_mr_cons_1 <- rma.mv(yi = lnRR,
                V = VCV_prey, 
                mods = ~ Log_area + Number_pattern + Treatment_stimulus,
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_prey)

summary(prey_mr_cons_1)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-29.8408   59.6816   75.6816   99.7131   76.7102   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0416  0.2039     16     no           Study_ID 
sigma^2.2  0.0142  0.1190     59     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    130     no          Cohort_ID 
sigma^2.4  0.0225  0.1498    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 149) = 605.0429, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 149) = 16.5329, p-val < .0001

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.3998  0.1980  -2.0193  149  0.0452  -0.7910 
Log_area                     0.2082  0.0471   4.4210  149  <.0001   0.1151 
Number_pattern              -0.0351  0.0177  -1.9773  149  0.0499  -0.0702 
Treatment_stimuluseyespot    0.1108  0.0643   1.7238  149  0.0868  -0.0162 
                             ci.ub      
intrcpt                    -0.0086    * 
Log_area                    0.3013  *** 
Number_pattern             -0.0000    * 
Treatment_stimuluseyespot   0.2379    . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(prey_mr_cons_1)

   R2_marginal R2_conditional 
     0.3411384      0.8107716 

Moderators associated with prey in preyator dataset

Code

prey_mr_pre <- rma.mv(yi = lnRR,
                V = VCV_prey, 
                mods = ~ Treatment_stimulus + Type_prey,
                random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_prey)

summary(prey_mr_pre)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-48.7681   97.5361  111.5361  132.6106  112.3249   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0836  0.2891     16     no           Study_ID 
sigma^2.2  0.0296  0.1719     59     no  Shared_control_ID 
sigma^2.3  0.0000  0.0001    130     no          Cohort_ID 
sigma^2.4  0.0293  0.1710    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 150) = 881.8004, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 150) = 1.6456, p-val = 0.1964

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                      0.2240  0.1109   2.0193  150  0.0452   0.0048 
Treatment_stimuluseyespot    0.1208  0.0750   1.6108  150  0.1093  -0.0274 
Type_preyreal               -0.2001  0.1771  -1.1294  150  0.2605  -0.5501 
                            ci.ub    
intrcpt                    0.4432  * 
Treatment_stimuluseyespot  0.2691    
Type_preyreal              0.1500    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(prey_mr_pre)

   R2_marginal R2_conditional 
    0.04918021     0.80464317 

3.4 Publication bias

Funnel plot

Code

# funnel plot - standard error
funnel(ma_prey, yaxis = "sei", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)" )

# funnel plot - inverse of standard error
prey_f2 <- funnel(ma_prey, yaxis = "seinv", 
      xlab = "Standarized residuals",
      ylab = "Precision (inverse of SE)",  col = c(alpha("orange", 0.5)))

Figure 23— Effect size and its standard error

Figure 24— Effect size and its inverse standard error

Figure 25— Funnel plot of effect size and its standard error

Egger’s test

Code

prey_df_bias <- dt_prey %>% mutate(sqrt_inv_e_n = sqrt((Cn + Tn)/(Cn * Tn)))

prey_bias_model <- rma.mv(yi = lnRR,
                      V = lnRR_var, 
                      mods = ~ 1 + sqrt_inv_e_n,
                      random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = prey_df_bias)

summary(prey_bias_model)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-50.6788  101.3576  113.3576  131.4613  113.9409   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0841  0.2900     16     no           Study_ID 
sigma^2.2  0.0248  0.1575     59     no  Shared_control_ID 
sigma^2.3  0.0027  0.0524    130     no          Cohort_ID 
sigma^2.4  0.0299  0.1728    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 896.6446, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 151) = 0.0310, p-val = 0.8604

Model Results:

              estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt         0.1978  0.1358  1.4564  151  0.1474  -0.0705  0.4662    
sqrt_inv_e_n    0.0784  0.4454  0.1761  151  0.8604  -0.8015  0.9584    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(prey_bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Figure 26— Egger’s test of publication bias

Decline effect

Code

prey_year_model <- rma.mv(yi = lnRR,
                      V = VCV_prey, 
                      mods = ~  1 + Year,
                      random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = prey_df_bias)

summary(prey_year_model)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-49.3996   98.7992  110.7992  128.9029  111.3825   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0599  0.2448     16     no           Study_ID 
sigma^2.2  0.0265  0.1628     59     no  Shared_control_ID 
sigma^2.3  0.0017  0.0418    130     no          Cohort_ID 
sigma^2.4  0.0306  0.1750    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 736.3783, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 151) = 2.8335, p-val = 0.0944

Model Results:

         estimate       se     tval   df    pval    ci.lb     ci.ub    
intrcpt   56.3493  33.3501   1.6896  151  0.0932  -9.5437  122.2423  . 
Year      -0.0279   0.0166  -1.6833  151  0.0944  -0.0607    0.0049  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(prey_year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

Figure 27— Decline effect

multi-moderator model

Code

prey_bias_all <- rma.mv(yi = lnRR,
                      V = VCV_prey, 
                      mods = ~ 1 + sqrt_inv_e_n + Year,
                      random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = prey_df_bias)

summary(prey_bias_all)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-48.8454   97.6909  111.6909  132.7653  112.4796   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0574  0.2396     16     no           Study_ID 
sigma^2.2  0.0283  0.1683     59     no  Shared_control_ID 
sigma^2.3  0.0013  0.0362    130     no          Cohort_ID 
sigma^2.4  0.0309  0.1757    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 150) = 730.1643, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 150) = 1.5150, p-val = 0.2232

Model Results:

              estimate       se     tval   df    pval    ci.lb     ci.ub    
intrcpt        58.6349  33.5654   1.7469  150  0.0827  -7.6872  124.9570  . 
sqrt_inv_e_n   -0.1440   0.4335  -0.3323  150  0.7401  -1.0006    0.7125    
Year           -0.0290   0.0167  -1.7405  150  0.0838  -0.0620    0.0039  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.5 Additional analyses

3.5.1 Background area (mm²)

Code

prey_mr_background <- rma.mv(yi = lnRR,
                        V = VCV_prey, 
                        mods = ~ Log_background,
                        random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                        test = "t",
                        method = "REML",
                        sparse = TRUE,
                        data = dt_prey)

summary(prey_mr_background)

Multivariate Meta-Analysis Model (k = 153; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-50.2053  100.4107  112.4107  130.5144  112.9940   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0779  0.2790     16     no           Study_ID 
sigma^2.2  0.0245  0.1564     59     no  Shared_control_ID 
sigma^2.3  0.0026  0.0509    130     no          Cohort_ID 
sigma^2.4  0.0304  0.1742    153     no             Obs_ID 

Test for Residual Heterogeneity:
QE(df = 151) = 905.4015, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 151) = 0.9913, p-val = 0.3210

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          -0.4590  0.6823  -0.6727  151  0.5022  -1.8069  0.8890    
Log_background    0.1004  0.1008   0.9956  151  0.3210  -0.0988  0.2996    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(prey_mr_background,
            mod = "Log_background",
            group = "Study_ID",
            xlab = "Log-transformed backgroud")

Figure 28— Effect size and background area

3.6 Figure galary

:::{.panel-tabset}

3.6.1 Discrete moderators

Code

# overall effect
prey_p1 <- orchard_plot(ma_prey,
                    group = "Study_ID",
                    xlab = "log response ratio (lnRR)", 
                    angle = 45,
                    legend.pos = "bottom.right") + 
                    scale_x_discrete(labels = c("Overall effect")) + 
                    theme(legend.direction = "horizontal",
                    legend.text = element_text(size = 6)) +
                    
                    scale_colour_okabe_ito()+
                    scale_fill_okabe_ito()

#eyespot vs. conspicuous patterns
prey_p2 <- orchard_plot(prey_mr_eyespot,
                  mod = "Treatment_stimulus",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6)) +
                  scale_colour_okabe_ito(order = 2:3)+
                  scale_fill_okabe_ito(order = 2:3)
                  
# prey type - real vs. artificial
prey_p3 <- orchard_plot(prey_mr_prey_type,
                  mod = "Type_prey",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6))+ 
                  scale_colour_okabe_ito(order = 4:5) +
                  scale_fill_okabe_ito(order = 4:5)

# combine
prey_patch1 <- prey_p1 | (prey_p2 / prey_p3)
prey_patch1 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig1.pdf", width = 10, height = 5, dpi = 450)

Figure 29— Discrete moderators

3.6.2 Continuous moderators

Code

# these figs were created by multi meta-regression model results
# log-transformed area
prey_p4 <- bubble_plot(prey_mr_area,
                  mod = "Log_area",
                  group = "Study_ID",
                  xlab = "Log-transformed area") +
                  scale_x_continuous(breaks = seq(0,7,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# number of patterns
prey_p5 <- bubble_plot(prey_mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") +
                  scale_x_continuous(breaks = seq(0,11,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# combine
prey_p4 / prey_p5 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig2_multi.pdf", width = 10, height = 10, dpi = 450)

Figure 30— Continuous moderators

3.6.3 Publication bias

Code

# funnel plot

# pdf("fig3.pdf")

funnel(ma_prey, yaxis = "seinv", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)",
      xlim = c(-4.0, 4.5), 
      ylim = c(0.01, 80.0), 
      col = c(alpha("mediumvioletred", 0.5)),
      cex = 0.7)

# dev.off()

# Egger's test and decline effect
prey_p7 <- bubble_plot(prey_bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size") +
            scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5)) 

prey_p8 <- bubble_plot(prey_year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication") +  scale_y_continuous(breaks = seq(-2.5, 4.0, 1.5)) 

# combine plots created by bubble_plot()
prey_pub <- prey_p7 / prey_p8
prey_pub +  plot_annotation(tag_levels = 'a') 

# output figure as pdf
# ggsave("fig4.pdf", width = 10, height = 10, dpi = 450)

(a) Funnel plot

(b) Egger’s test and Decline effect

Figure 31— Publication bias

4 Meta-analysis - all dataset

We used meta-analytical models to calculate total I2 (a measure of heterogeneity not caused by sampling error; Higgins. et al., 2003) and the partial I2 explained by each random factor.

Which random effects remove in our models?

We can remove Shared control ID and Cohort ID

Code

# I exclude cohort_ID because sigma^2.2 = 0 and I2 = 0
ma_all <- rma.mv(yi = lnRR,
                  V = VCV, 
                  random = list(~1 | Study_ID,
                                ~1 | Shared_control_ID,
                                ~1 | Cohort_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(ma_all)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-251.0417   502.0833   512.0833   530.0569   512.3115   

Variance Components:

            estim    sqrt  nlvls  fixed             factor 
sigma^2.1  0.0485  0.2203     33     no           Study_ID 
sigma^2.2  0.0000  0.0000     89     no  Shared_control_ID 
sigma^2.3  0.0000  0.0000    164     no          Cohort_ID 
sigma^2.4  0.2334  0.4831    270     no             Obs_ID 

Test for Heterogeneity:
Q(df = 269) = 2729.8913, p-val < .0001

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.1977  0.0581  3.4039  269  0.0008  0.0833  0.3120  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

i2 <- round(i2_ml(ma_all), 4)
r2 <- round(r2_ml(ma_all), 4)

I2_Total	I2_Study_ID	I2_Shared_control_ID	I2_Cohort_ID	I2_Obs_ID
96.4974	16.6131	0	0	79.8844

R2_marginal	R2_conditional
0	0.1722

Table 3 - Heterogeneity of effect size

Code

orchard_plot(ma_all,
              group = "Study_ID",
              xlab = "log response ratio (lnRR)", 
              angle = 45) + 
              scale_x_discrete(labels = c("Overall effect")) + 
              scale_colour_okabe_ito()+
              scale_fill_okabe_ito()

Figure 32— Estimates of overall effect size and 95% confidence intervals

4.1 Meta-regressions: uni-moderator

Pattern type

Eyespot vs. conspicuous patterns - Is there significant difference of effect size between two patterns?

Code

#normal model
mr_eyespot <- rma.mv(yi = lnRR,
                    V = VCV, 
                    mods = ~ Treatment_stimulus,
                    random = list(~1 | Study_ID,
                                  ~1 | Obs_ID),
                    test = "t",
                    method = "REML", 
                    sparse = TRUE,
                    data = dt_all)

summary(mr_eyespot)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-250.1918   500.3835   508.3835   522.7475   508.5356   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0526  0.2293     33     no  Study_ID 
sigma^2.2  0.2335  0.4833    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2721.6080, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.3291, p-val = 0.5666

Model Results:

                           estimate      se    tval   df    pval    ci.lb 
intrcpt                      0.1579  0.0921  1.7132  268  0.0878  -0.0236 
Treatment_stimuluseyespot    0.0603  0.1050  0.5737  268  0.5666  -0.1465 
                            ci.ub    
intrcpt                    0.3393  . 
Treatment_stimuluseyespot  0.2671    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_1 <- round(r2_ml(mr_eyespot), 4)

R2_marginal	R2_conditional
0.0027	0.186

Table XX. Model goodness-of-fit

Code

orchard_plot(mr_eyespot,
            mod = "Treatment_stimulus",
            group = "Study_ID",
            xlab = "log response ratio (lnRR)",
            angle = 45) + 
            scale_colour_okabe_ito(order = 2:3)+
            scale_fill_okabe_ito(order = 2:3)

Figure 33— Eyespot vs. conspicuous patterns

Code

# intercept-removed model
mr_eyespot1 <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~ Treatment_stimulus -1,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_all)

summary(mr_eyespot1)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-250.1918   500.3835   508.3835   522.7475   508.5356   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0526  0.2293     33     no  Study_ID 
sigma^2.2  0.2335  0.4833    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2721.6080, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 268) = 5.7321, p-val = 0.0037

Model Results:

                               estimate      se    tval   df    pval    ci.lb 
Treatment_stimulusconspicuous    0.1579  0.0921  1.7132  268  0.0878  -0.0236 
Treatment_stimuluseyespot        0.2181  0.0688  3.1722  268  0.0017   0.0827 
                                ci.ub     
Treatment_stimulusconspicuous  0.3393   . 
Treatment_stimuluseyespot      0.3535  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pattern diameter (mm)

Code

mr_diameter <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~ Log_diameter,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = dt_all)

summary(mr_diameter)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-247.2115   494.4230   502.4230   516.7869   502.5751   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0440  0.2099     33     no  Study_ID 
sigma^2.2  0.2286  0.4781    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2720.8294, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 6.0336, p-val = 0.0147

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt        -0.1757  0.1618  -1.0858  268  0.2785  -0.4942  0.1428    
Log_diameter    0.1931  0.0786   2.4563  268  0.0147   0.0383  0.3479  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_2 <- round(r2_ml(mr_diameter), 4)

R2_marginal	R2_conditional
0.0662	0.2171

Table XX. Model goodness-of-fit

Code

bubble_plot(mr_diameter,
            mod = "Log_diameter",
            group = "Study_ID",
            xlab = "Log-transformed diameter")

Pattern area (mm²)

Code

mr_area <- rma.mv(yi = lnRR,
                  V = VCV, 
                  mods = ~ Log_area,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(mr_area)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-246.5105   493.0210   501.0210   515.3849   501.1731   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0472  0.2172     33     no  Study_ID 
sigma^2.2  0.2265  0.4759    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2718.2726, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 7.3561, p-val = 0.0071

Model Results:

          estimate      se     tval   df    pval    ci.lb   ci.ub     
intrcpt    -0.1972  0.1563  -1.2616  268  0.2082  -0.5049  0.1105     
Log_area    0.1119  0.0412   2.7122  268  0.0071   0.0307  0.1931  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_3 <- round(r2_ml(mr_area), 4)

R2_marginal	R2_conditional
0.0856	0.2433

Table XX. Model goodness-of-fit

Code

bubble_plot(mr_area,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Pattern number

Code

mr_num <- rma.mv(yi = lnRR,
                  V = VCV,
                  mods = ~ Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

summary(mr_num)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-248.0452   496.0904   504.0904   518.4544   504.2425   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0442  0.2102     33     no  Study_ID 
sigma^2.2  0.2296  0.4792    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2703.4427, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 5.2375, p-val = 0.0229

Model Results:

                estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt           0.3437  0.0854   4.0229  268  <.0001   0.1755   0.5119  *** 
Number_pattern   -0.0584  0.0255  -2.2886  268  0.0229  -0.1086  -0.0082    * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_4 <- round(r2_ml(mr_num), 4)

R2_marginal	R2_conditional
0.0246	0.182

Table XX. Model goodness-of-fit

Code

bubble_plot(mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") 

Background area (mm²)

Code

mr_background <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Log_background,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                        test = "t",
                        method = "REML",
                        sparse = TRUE,
                        data = dt_all)

summary(mr_background)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-249.8534   499.7068   507.7068   522.0707   507.8589   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0545  0.2334     33     no  Study_ID 
sigma^2.2  0.2331  0.4828    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2667.6510, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.2291, p-val = 0.6326

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt           0.3997  0.4246   0.9413  268  0.3474  -0.4363  1.2356    
Log_background   -0.0289  0.0604  -0.4786  268  0.6326  -0.1479  0.0900    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_background)

   R2_marginal R2_conditional 
    0.00417208     0.19276397 

Code

bubble_plot(mr_background,
            mod = "Log_background",
            group = "Study_ID",
            xlab = "Log-transformed backgroud")

Figure 34— Effect size and background area

Prey material type

Our dataset includes two material types of prey: live and artificial. Is there significant difference of effect size between two stimuli?

Code

# normal model
mr_prey_type <- rma.mv(yi = lnRR,
                  V = VCV, 
                  mods = ~ Type_prey,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML", 
                  sparse = TRUE,
                  data = dt_all)

summary(mr_prey_type)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-250.2013   500.4026   508.4026   522.7666   508.5547   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0552  0.2351     33     no  Study_ID 
sigma^2.2  0.2329  0.4826    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2702.2195, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.1246, p-val = 0.7244

Model Results:

               estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt          0.1829  0.0749  2.4401  268  0.0153   0.0353  0.3304  * 
Type_preyreal    0.0444  0.1259  0.3529  268  0.7244  -0.2035  0.2924    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_5 <- round(r2_ml(mr_prey_type), 4)

R2_marginal	R2_conditional
0.0013	0.1928

Table XX. Model goodness-of-fit

Code

orchard_plot(mr_prey_type,
              mod = "Type_prey",
              group = "Study_ID",
              xlab = "Type of prey",
              angle = 45) +
            scale_colour_okabe_ito(order = 4:5)+
            scale_fill_okabe_ito(order = 4:5)

Code

# intercept-removed model
mr_prey_type1 <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Type_prey -1,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                        test = "t",
                        method = "REML", 
                        sparse = TRUE,
                        data = dt_all)

summary(mr_prey_type1)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-250.2013   500.4026   508.4026   522.7666   508.5547   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0552  0.2351     33     no  Study_ID 
sigma^2.2  0.2329  0.4826    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2702.2195, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 268) = 5.4986, p-val = 0.0046

Model Results:

                     estimate      se    tval   df    pval   ci.lb   ci.ub    
Type_preyartificial    0.1829  0.0749  2.4401  268  0.0153  0.0353  0.3304  * 
Type_preyreal          0.2273  0.1012  2.2457  268  0.0255  0.0280  0.4266  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prey shape type

Our dataset includes 4 types of prey type: real butterfly, abstract butterfly, abstract caterpillar, and abstract prey. Is there significant difference of effect size between two stimuli?

Code

mr_prey_shape <- rma.mv(yi = lnRR,
                        V = VCV, 
                        mods = ~ Shape_prey -1,
                        random = list(~1 | Study_ID,
                                      ~1 | Obs_ID),
                                      test = "t",
                                      method = "REML", 
                                      sparse = TRUE,
                                      data = dt_all)

summary(mr_prey_shape)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-246.3898   492.7797   504.7797   526.2806   505.1040   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0537  0.2317     33     no  Study_ID 
sigma^2.2  0.2305  0.4801    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 266) = 2509.2480, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 266) = 3.6583, p-val = 0.0064

Model Results:

                                estimate      se    tval   df    pval    ci.lb 
Shape_preyabstract_butterfly      0.3219  0.1057  3.0446  266  0.0026   0.1137 
Shape_preyabstract_caterpillar    0.0661  0.1259  0.5246  266  0.6003  -0.1819 
Shape_preyabstract_prey           0.0113  0.1842  0.0615  266  0.9510  -0.3514 
Shape_preybutterfly               0.2263  0.1004  2.2549  266  0.0250   0.0287 
                                 ci.ub     
Shape_preyabstract_butterfly    0.5301  ** 
Shape_preyabstract_caterpillar  0.3140     
Shape_preyabstract_prey         0.3741     
Shape_preybutterfly             0.4239   * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_6 <- round(r2_ml(mr_prey_shape), 4)

dat_all$Shape_prey <- factor(dat_all$Shape_prey)
levels(dat_all$Shape_prey)

[1] "abstract_butterfly"   "abstract_caterpillar" "abstract_prey"       
[4] "butterfly"           

Code

mat_ex <- 
    cbind(contrMat(rep(1,
                       length(mr_prey_shape$ci.ub)),
                   type = "Tukey"))
sig_test <- summary(glht(mr_prey_shape, linfct= mat_ex), test = adjusted("none"))

sig_test

     Simultaneous Tests for General Linear Hypotheses

Fit: rma.mv(yi = lnRR, V = VCV, mods = ~Shape_prey - 1, random = list(~1 | 
    Study_ID, ~1 | Obs_ID), data = dt_all, method = "REML", test = "t", 
    sparse = TRUE)

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0 -0.25583    0.16444  -1.556    0.120
3 - 1 == 0 -0.31057    0.21242  -1.462    0.144
4 - 1 == 0 -0.09561    0.14577  -0.656    0.512
3 - 2 == 0 -0.05474    0.22317  -0.245    0.806
4 - 2 == 0  0.16022    0.16104   0.995    0.320
4 - 3 == 0  0.21496    0.20980   1.025    0.306
(Adjusted p values reported -- none method)

R2_marginal	R2_conditional
0.0521	0.2312

Table XX. Model goodness-of-fit

Code

orchard_plot(mr_prey_shape,
              mod = "Shape_prey",
              group = "Study_ID",
              xlab = "Shape of prey",
              angle = 45) +
              scale_colour_okabe_ito(order = 6:9)+
              scale_fill_okabe_ito(order = 6:9)

Figure 35— Effect size and prey shape

4.2 Correlation visualisation and choose moderators

Before we run multi-moderator meta-regressions, we need to consider the correlation between moderators. Area, diameter and background seem to be correlated. We visualised the correlation between these variables.

Visualisation of correlation

Code

# correlation between continuous variables
corr_cont <- round(cor(dt_all[, c("Diameter_pattern", "Area_pattern", "Number_pattern", "Area_background")]),2)

p1 <- ggcorrplot(corr_cont, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(a) Continuous variables")

corr_cont_log <- round(cor(dt_all[, c("Log_diameter", "Log_area", "Number_pattern", "Log_background")]),2)

p2 <- ggcorrplot(corr_cont_log, hc.order = TRUE, lab = TRUE, outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"), title = "(b) Log-transormed continuous variables")

p1 + p2 +plot_layout(guides = 'collect')

Figure 36— Correlation between coninuous moderators

Code

dat1 <- dt_all %>%
  dplyr::select("Treatment_stimulus", "Type_prey", "Shape_prey")

corr_cat <- GKtauDataframe(dat1)
plot(corr_cat)

Figure 37— Correlation between categorical variables

We should not include “Shape_prey” and “Type_prey” in the model at the same time. Because they are correlated.

Choose moderators

We used model R2 values to find better model and modelator VIF values to check multicollinearity between moderators. Higher R2 indicates better model and VIF > 2 indicates multicollinearity.

4.2.0.1 R2

Code

r2_area <- rma.mv(yi = lnRR,
                  V = VCV,
                  mods = ~ Log_area + Number_pattern,
                  random = list(~1 | Study_ID,
                                ~1 | Obs_ID),
                  test = "t",
                  method = "REML",
                  sparse = TRUE,
                  data = dt_all)

r2_diameter <- rma.mv(yi = lnRR,
                      V = VCV,
                      mods = ~ Log_diameter + Number_pattern,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

# check r2 values
r2_ml(r2_area)

   R2_marginal R2_conditional 
    0.08098982     0.22938197 

Code

r2_ml(r2_diameter)

   R2_marginal R2_conditional 
    0.06846176     0.20934212 

It seems area is a better predictor than diameter .

4.3 Meta-regressions: multi-moderator

All potential moderators (full model)

Include all moderators in the model.

Code

mr_all <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Treatment_stimulus + Log_area 
                + Number_pattern + Type_prey,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_all)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.6478   485.2957   499.2957   524.3538   499.7315   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0522  0.2285     33     no  Study_ID 
sigma^2.2  0.2244  0.4738    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 265) = 2598.7375, p-val < .0001

Test of Moderators (coefficients 2:5):
F(df1 = 4, df2 = 265) = 2.8780, p-val = 0.0233

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.0851  0.2100  -0.4053  265  0.6856  -0.4987 
Treatment_stimuluseyespot    0.0216  0.1065   0.2027  265  0.8395  -0.1881 
Log_area                     0.0964  0.0446   2.1612  265  0.0316   0.0086 
Number_pattern              -0.0527  0.0287  -1.8393  265  0.0670  -0.1092 
Type_preyreal                0.1776  0.1358   1.3080  265  0.1920  -0.0898 
                            ci.ub    
intrcpt                    0.3284    
Treatment_stimuluseyespot  0.2312    
Log_area                   0.1841  * 
Number_pattern             0.0037  . 
Type_preyreal              0.4450    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_all)

   R2_marginal R2_conditional 
    0.08328291     0.25633953 

Moderators associated with conspicuousness

Only include moderators related to conspicuousness (pattern area and number of patterns) in the model.

Code

mr_cons <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Log_area + Number_pattern,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_cons)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-244.9501   489.9002   499.9002   517.8365   500.1301   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0436  0.2087     33     no  Study_ID 
sigma^2.2  0.2262  0.4756    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 267) = 2679.2370, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 267) = 4.8661, p-val = 0.0084

Model Results:

                estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt          -0.0185  0.1916  -0.0966  267  0.9231  -0.3959  0.3588    
Log_area          0.0899  0.0427   2.1056  267  0.0362   0.0058  0.1740  * 
Number_pattern   -0.0406  0.0267  -1.5215  267  0.1293  -0.0932  0.0119    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_cons)

   R2_marginal R2_conditional 
    0.08098982     0.22938197 

Code

bubble_plot(mr_cons,
            mod = "Log_area",
            group = "Study_ID",
            xlab = "Log-transformed area")

Code

mr_cons_1 <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Log_area + Number_pattern + Treatment_stimulus,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_cons_1)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-244.1838   488.3676   500.3676   521.8686   500.6920   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0480  0.2190     33     no  Study_ID 
sigma^2.2  0.2262  0.4756    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 266) = 2679.0362, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 266) = 3.2831, p-val = 0.0214

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
intrcpt                     -0.0540  0.2057  -0.2625  266  0.7931  -0.4590 
Log_area                     0.0907  0.0437   2.0761  266  0.0388   0.0047 
Number_pattern              -0.0405  0.0269  -1.5039  266  0.1338  -0.0935 
Treatment_stimuluseyespot    0.0506  0.1025   0.4939  266  0.6218  -0.1512 
                            ci.ub    
intrcpt                    0.3510    
Log_area                   0.1767  * 
Number_pattern             0.0125    
Treatment_stimuluseyespot  0.2524    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_cons_1)

   R2_marginal R2_conditional 
    0.08567312     0.24568295 

Moderators associated with prey

Code

mr_pre <- rma.mv(yi = lnRR,
                V = VCV, 
                mods = ~ Treatment_stimulus + Type_prey,
                random = list(~1 | Study_ID,
                              ~1 | Obs_ID),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dt_all)

summary(mr_pre)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-249.3570   498.7140   508.7140   526.6503   508.9439   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0589  0.2426     33     no  Study_ID 
sigma^2.2  0.2333  0.4830    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 267) = 2690.3970, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 267) = 0.2046, p-val = 0.8151

Model Results:

                           estimate      se    tval   df    pval    ci.lb 
intrcpt                      0.1487  0.0998  1.4907  267  0.1372  -0.0477 
Treatment_stimuluseyespot    0.0569  0.1093  0.5204  267  0.6032  -0.1584 
Type_preyreal                0.0338  0.1308  0.2580  267  0.7966  -0.2238 
                            ci.ub    
intrcpt                    0.3452    
Treatment_stimuluseyespot  0.2721    
Type_preyreal              0.2913    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

r2_ml(mr_pre)

   R2_marginal R2_conditional 
   0.003607887    0.204348070 

4.4 Publication bias

Funnel plot

Code

# funnel plot - standard error
funnel(ma_all, yaxis = "sei", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)" )

# funnel plot - inverse of standard error
f2 <- funnel(ma_all, yaxis = "seinv", 
      xlab = "Standarized residuals",
      ylab = "Precision (inverse of SE)",  col = c(alpha("orange", 0.5)))

Figure 38— Effect size and its standard error

Figure 39— Effect size and its inverse standard error

Figure 40— Funnel plot of effect size and its standard error

Egger’s test

Code

df_bias <- dt_all %>% mutate(sqrt_inv_e_n = sqrt((Cn + Tn)/(Cn * Tn)))

bias_model <- rma.mv(yi = lnRR,
                      V = lnRR_var, 
                      mods = ~ 1 + sqrt_inv_e_n,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(bias_model)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-248.4547   496.9094   504.9094   519.2734   505.0615   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0914  0.3024     33     no  Study_ID 
sigma^2.2  0.2074  0.4554    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2155.7534, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.0568, p-val = 0.8118

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt         0.2491  0.1353   1.8404  268  0.0668  -0.0174  0.5155  . 
sqrt_inv_e_n   -0.0895  0.3757  -0.2383  268  0.8118  -0.8291  0.6501    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Figure 41— Egger’s test of publication bias

Decline effect

Code

year_model <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~  1 + Year,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(year_model)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-250.1678   500.3355   508.3355   522.6995   508.4876   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0510  0.2259     33     no  Study_ID 
sigma^2.2  0.2339  0.4836    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2659.1976, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.0147, p-val = 0.9035

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt    1.7561  12.8324   0.1368  268  0.8913  -23.5091  27.0212    
Year      -0.0008   0.0064  -0.1214  268  0.9035   -0.0134   0.0118    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

Figure 42— Decline effect

Multi-moderator

Code

multi_bias <- rma.mv(yi = lnRR,
                      V = VCV, 
                      mods = ~sqrt_inv_e_n + Year,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML", 
                      sparse = TRUE,
                      data = df_bias)

summary(multi_bias)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-248.8289   497.6578   507.6578   525.5940   507.8876   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0415  0.2038     33     no  Study_ID 
sigma^2.2  0.2360  0.4858    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 267) = 2599.6744, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 267) = 0.5923, p-val = 0.5538

Model Results:

              estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt         5.7530  12.7722   0.4504  267  0.6528  -19.3941  30.9000    
sqrt_inv_e_n   -0.3953   0.3686  -1.0724  267  0.2845   -1.1211   0.3304    
Year           -0.0027   0.0063  -0.4266  267  0.6700   -0.0152   0.0098    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(multi_bias,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size")

Code

bubble_plot(multi_bias,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication")

4.5 Additional analyses

4.5.1 Dataset

We compared effect size between two datasets (predator and prey) to see whether there is any difference in effect size between two datasets.

Code

mr_dataset <- rma.mv(yi = lnRR,
                      V = VCV,
                      mods = ~ Dataset,
                      random = list(~1 | Study_ID,
                                    ~1 | Obs_ID),
                      test = "t",
                      method = "REML",
                      sparse = TRUE,
                      data = dt_all)

summary(mr_dataset)

Multivariate Meta-Analysis Model (k = 270; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-249.9719   499.9439   507.9439   522.3078   508.0959   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0519  0.2279     33     no  Study_ID 
sigma^2.2  0.2331  0.4828    270     no    Obs_ID 

Test for Residual Heterogeneity:
QE(df = 268) = 2658.6695, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 268) = 0.4221, p-val = 0.5164

Model Results:

             estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt        0.1563  0.0875  1.7852  268  0.0754  -0.0161  0.3286  . 
Datasetprey    0.0756  0.1164  0.6497  268  0.5164  -0.1535  0.3047    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mr_dataset,
            mod = "Dataset",
            group = "Study_ID",
            xlab = "Dataset", 
            angle = 45) 

Figure 43— Effect size and dataset - predator and prey

4.6 Figure galary

Discrete moderators

Code

# overall effect
p1 <- orchard_plot(ma_all,
                    group = "Study_ID",
                    xlab = "log response ratio (lnRR)", 
                    angle = 45,
                    legend.pos = "bottom.right") + 
                    scale_x_discrete(labels = c("Overall effect")) + 
                    theme(legend.direction = "horizontal",
                    legend.text = element_text(size = 6)) +
                    
                    scale_colour_okabe_ito()+
                    scale_fill_okabe_ito()

#eyespot vs. conspicuous patterns
p2 <- orchard_plot(mr_eyespot,
                  mod = "Treatment_stimulus",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6)) +
                  scale_colour_okabe_ito(order = 2:3)+
                  scale_fill_okabe_ito(order = 2:3)
                  
# prey type - real vs. artificial
p3 <- orchard_plot(mr_prey_type,
                  mod = "Type_prey",
                  group = "Study_ID",
                  xlab = "log response ratio (lnRR)",
                  angle = 45,
                  legend.pos = "bottom.right") +
                  theme(legend.direction = "horizontal",
                  legend.text = element_text(size = 6))+ 
                  scale_colour_okabe_ito(order = 4:5) +
                  scale_fill_okabe_ito(order = 4:5)

# combine
patch1 <- p1 | (p2 / p3)
patch1 + plot_annotation(tag_levels = "a")

# output figure as pdf
# ggsave("fig1.pdf", width = 10, height = 5, dpi = 450)

Figure 44— Discrete moderators

Continuous moderators

Code

# these figs were created by multi meta-regression model results
# log-transformed area
p4 <- bubble_plot(mr_area,
                  mod = "Log_area",
                  group = "Study_ID",
                  xlab = "Log-transformed area") +
                  scale_x_continuous(breaks = seq(0,7,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# number of patterns
p5 <- bubble_plot(mr_num,
            mod = "Number_pattern",
            group = "Study_ID",
            xlab = "Number of patterns") +
                  scale_x_continuous(breaks = seq(0,11,1.5)) +
                  scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5))

# combine
p4 / p5 + plot_annotation(tag_levels = "a")

# output figure as pdf
ggsave("fig2_multi.pdf", dpi = 450, units = "mm")

Figure 45— Continuous moderators

Publication bias

Code

# funnel plot

# pdf("fig3.pdf")

funnel(ma_all, yaxis = "seinv", 
      xlab = "Standarised residuals",
      ylab = "Precision (inverse of SE)",
      xlim = c(-4.0, 4.5), 
      ylim = c(0.01, 80.0), 
      col = c(alpha("mediumvioletred", 0.5)),
      cex = 0.7)

# dev.off()

# Egger's test and decline effect
p7 <- bubble_plot(bias_model,
            mod = "sqrt_inv_e_n",
            group = "Study_ID",
            xlab = "Square root of inverse of effective sample size") +
            scale_y_continuous(breaks = seq(-2.5, 4.5, 1.5)) 

p8 <- bubble_plot(year_model,
            mod = "Year",
            group = "Study_ID",
            xlab = "Year of publication") +  scale_y_continuous(breaks = seq(-2.5, 4.0, 1.5)) 

# combine plots created by bubble_plot()
pub <- p7 / p8
pub +  plot_annotation(tag_levels = 'a') 

# output figure as pdf
# ggsave("fig4.pdf", width = 10, height = 10, dpi = 450)

Figure 46— Funnel plot

Figure 47— Egger’s test and Decline effect

Figure 48— Publication bias

5 R session infomation

R version 4.3.1 (2023-06-16)

Platform: aarch64-apple-darwin20 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: parallel, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: orchaRd(v.2.0), multcomp(v.1.4-25), TH.data(v.1.1-2), MASS(v.7.3-60), survival(v.3.5-5), mvtnorm(v.1.2-3), miWQS(v.0.4.4), metafor(v.4.2-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.6-1), patchwork(v.1.1.3), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.3), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), tidyverse(v.2.0.0), knitr(v.1.43), here(v.1.0.1), ggcorrplot(v.0.1.4), ggplot2(v.3.4.3), ggokabeito(v.0.1.0), dtplyr(v.1.3.1), DT(v.0.29), GoodmanKruskal(v.0.0.3) and ape(v.5.7-1)

loaded via a namespace (and not attached): gridExtra(v.2.3), glm2(v.1.2.1), sandwich(v.3.0-2), rlang(v.1.1.1), magrittr(v.2.0.3), compiler(v.4.3.1), mgcv(v.1.8-42), systemfonts(v.1.0.4), reshape2(v.1.4.4), vctrs(v.0.6.3), quantreg(v.5.97), pkgconfig(v.2.0.3), fastmap(v.1.1.1), ellipsis(v.0.3.2), backports(v.1.4.1), mcmc(v.0.9-7), labeling(v.0.4.3), pander(v.0.6.5), utf8(v.1.2.3), rmarkdown(v.2.24), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), ragg(v.1.2.5), MatrixModels(v.0.5-2), xfun(v.0.40), cachem(v.1.0.8), jsonlite(v.1.8.7), gmm(v.1.8), cluster(v.2.1.4), R6(v.2.5.1), bslib(v.0.5.1), stringi(v.1.7.12), Rsolnp(v.1.16), rlist(v.0.4.6.2), rpart(v.4.1.19), estimability(v.1.4.1), jquerylib(v.0.1.4), Rcpp(v.1.0.11), zoo(v.1.8-12), tmvmixnorm(v.1.1.1), base64enc(v.0.1-3), pacman(v.0.5.1), splines(v.4.3.1), nnet(v.7.3-19), timechange(v.0.2.0), tidyselect(v.1.2.0), rstudioapi(v.0.15.0), yaml(v.2.3.7), codetools(v.0.2-19), plyr(v.1.8.8), lattice(v.0.21-8), withr(v.2.5.0), coda(v.0.19-4), evaluate(v.0.21), tmvtnorm(v.1.5), foreign(v.0.8-84), pillar(v.1.9.0), corrplot(v.0.92), checkmate(v.2.2.0), stats4(v.4.3.1), generics(v.0.1.3), rprojroot(v.2.0.3), invgamma(v.1.1), mathjaxr(v.1.6-0), truncnorm(v.1.0-9), hms(v.1.1.3), munsell(v.0.5.0), scales(v.1.2.1), xtable(v.1.8-4), glue(v.1.6.2), emmeans(v.1.8.8), Hmisc(v.5.1-0), tools(v.4.3.1), data.table(v.1.14.8), SparseM(v.1.81), matrixNormal(v.0.1.1), grid(v.4.3.1), MCMCpack(v.1.6-3), crosstalk(v.1.2.0), colorspace(v.2.1-0), nlme(v.3.1-162), beeswarm(v.0.4.0), condMVNorm(v.2020.1), vipor(v.0.4.5), htmlTable(v.2.4.1), latex2exp(v.0.9.6), Formula(v.1.2-5), cli(v.3.6.1), textshaping(v.0.3.6), fansi(v.1.0.4), gtable(v.0.3.4), sass(v.0.4.7), digest(v.0.6.33), farver(v.2.1.1), htmlwidgets(v.1.6.2), htmltools(v.0.5.6) and lifecycle(v.1.0.3)